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Chapter I 


Introduction and summary 


Seismological research has provided a major contribution to our present knowledge of 
the structure of the Earth. In the 1930’s, the radially symmetric seismic velocity 
distribution was established, showing regions of different seismic character (crust, upper 
mantle, transition zone, lower mantle, inner and outer core). The interpretation of the 
seismological findings in terms of the compositional, thermal, and geodynamical 
constitution of the Earth has, since then, been an important subject in solid Earth 
geophysics. Presently, there are still major outstanding questions in this field that require 
high resolution seismological data for their solution. The research of this thesis must be 
seen in this context: it addresses two aspects that concern the detailed structure of the 
upper mantle. One is related to the fine structure of the sharp seismic velocity increases 
( discontinuities’) at depths of about 400 and 670 km, marking the upper and lower 
boundary of the transition zone. The other concerns the degree of lateral (seismic velocity) 
heterogeneity in the upper mantle. In the following, a brief outline is presented of the 
discussion about the nature of the upper mantle discontinuities, and of several aspects 
associated with upper mantle lateral heterogeneity. 


A change in gradient of the first arrival travel time curve near 20° led to the discovery 
of the 20°, or 400-km discontinuity. After the first seismological observations of this 20° 
discontinuity, the transition was explained by a phase transformation (see Nolet & Wortel, 
1988, for a brief review). Although chemical changes in the transition region have been 
proposed more often (e.g. Anderson, 1968; Press, 1968), the first more detailed 
interpretation of the 400-km discontinuity as being of chemical nature was made by Bass & 
Anderson (1984) and Anderson & Bass (1984). They suggested that this boundary marks 
the transition from olivine-rich material (pyrolite) to more eclogite-rich material (piclogite). 
This latter composition would better fit the seismic velocities in the transition zone and the 
sharpness of the 400-km discontinuity, but their arguments have been disputed by others. 
Weidner (1985) and Bina & Wood (1987) showed that the olivine-B-spinel phase 
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transformation in a pyrolitic composition could explain the seismic observations, and that a 
chemical change at a depth of about 400 km need not be invoked. 

First evidence for the presence of the 670-km discontinuity was obtained in the 1960’s 
(e.g. Niazi & Anderson, 1965; Johnson, 1967; Green & Hales, 1968). Since then, seismic 
observations of refracted and reflected phases at this boundary have increased substantially, 
but, as with the 400-km discontinuity, there is still a debate about its nature. Experimental 
studies have shown that phase transformations in systems of the two most abundant 
constituents of the mantle, olivine and pyroxene, occur at pressure and temperature 
conditions similar to those in the Earth at a depth of about 650 km (e.g. Liu, 1977; Yagi et 
al., 1979). These results seem to point to a phase transition interpretation. However, the 
sharpness of this discontinuity, as suggested by seismological studies, would favour the 
hypothesis of a chemical boundary. 

As is clear from the foregoing, there is presently no general agreement about the 
nature of the two upper mantle discontinuities. 


Since the late 1960’s, it is well established that there are regional variations in the 
upper mantle structure. Surface wave studies, in particular, have provided a significant 
contribution to this knowledge, as Knopoff (1972) showed that a tectonic regionalization 
could be made on the basis of the dispersion of surface waves. In more recent surface 
wave studies, advanced processing techniques are employed to optimize the seismic 
velocity information acquired from the data. The NARS array in western Europe (chapter 
2) has provided suitable data for several of such investigations. Using higher mode surface 
waves, Dost (1987) determined the shear wave velocity and density distribution of the 
western European Platform and Scandinavian Shield, and showed that these regions have 
different shear velocities to a depth of 400 km. Nolet & Scheffers (1988) used nonlinear 
waveform fitting of higher modes surface waves to find deep lateral variations under the 
NARS array itself. Snieder (1987) used scattered and direct surface waves to derive a 
three-dimensional model of the shear-wave velocity structure to a depth of 200 km beneath 
Europe. 

Lateral variations in the upper mantle structure have also been inferred from body 
waves. After the first few regional body wave studies it was recognized that the upper 
mantle seismic velocity distribution varies from region to region (e.g. Archambeau et al., 
1969). These regional models present relatively accurate velocity-depth profiles, but do not 
yield information about the horizontal scale over which the variations occur. Tomography 
provides a technique to map these lateral variations, and tomographic studies of P delay 
times, have resulted in three-dimensional representations of the upper mantle (for Europe 
e.g.: Babuska et al. 1984; Spakman et al., 1988). 

All of these results lack, however, the resolution to verify different hypotheses 
proposed for the structure and evolution of the upper mantle. Ringwood & Irifune (1988) 
consider the major element composition of the mantle to be uniform on a large scale (>100 
km), while differentiation at mid-ocean ridges and subduction at trenches cause medium- 


Introduction and summary Il 


scale heterogeneity (5-100 km). Jordan (1975) assumes that the seismic velocity variations 
in the upper mantle are essentially correlated with age. Vlaar (1983) considers the 
implications of subhorizontal subduction of young oceanic lithosphere, which will have its 
imprint on the local upper mantle structure with expected variations on a length scale of 
100 km or less. 


It is clear that detailed seismological evidence may contribute to an improved insight 
into the constitution of the Earth’s upper mantle. Chapter 4 and 5 of this thesis present the 
results of two detailed body wave studies of the upper mantle that address the question of 
the degree of lateral heterogeneity and the fine structure of the upper mantle discontinuities, 
respectively. The data used for these studies are for the major part from stations of the 
NARS array in western Europe. A brief description of this seismic array is given in chapter 
2. Another introductory chapter (chapter 3) discusses the merits and shortcomings of the 
WKBJ-technique (Chapman, 1978): a method to calculate synthetic seismograms, that has 
been employed extensively in this thesis. 

Chapter 4 describes how the seismic velocity distribution of the upper mantle beneath 
Europe is investigated by body wave modelling of seismograms from earthquakes in the 
eastern Mediterranean recorded by stations of the NARS array. It is shown that the gross 
seismic velocity structure of the subcontinental European mantle can be described by a 
high velocity lid, a low velocity zone, and discontinuities at depths of about 420 and 650 
km. The data, however, also present evidence for distinct small scale lateral heterogeneity 
to a depth of about 400 km, indicating that the upper mantle of this area may be less 
homogeneous than was generally assumed. 

Evidence for a locally sharp 670-km discontinuity is inferred from waves that have 
undergone wavetype conversion at the upper mantle discontinuities. Chapter 5 presents the 
results of this study of P-to-S converted waves, for which not only NARS data have been 
used, but also short-period seismograms from the RSTN network in the United States. The 
results of this study may imply that the upper and lower mantle are chemically distinct. 

In the last chapter of this thesis (chapter 6) an overview is given of our insights into 
the upper mantle structure. The implications of the seismological results are presented in 
the light of evidence from other fields and of hypotheses proposed for the upper mantle. 
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Chapter 2 


The NARS array 


The Network of Autonomously Recording Stations (NARS) was primarily developed 
for dispersion measurements of higher mode surface waves across the western European 
Platform. However, for the NARS experiment in western Europe (1983-1987), the 
instruments were not only tuned for the recording of teleseismic surface waves (Dost, 
1987), but also for an adequate registration of body waves (data on which this thesis is 
based). 


2.1 NARS SYSTEM RESPONSE 


The system design of the NARS stations is displayed in fig. 2.1. The output of a 
vertical long-period Teledyne Geotech (SL 210) seismometer and two horizontal (SL 220) 
seismometers undergoes pre-amplification and electronic shaping before it is digitized in 
the Kinemetrics PDR-2 digital event recorder. The internal clock of the PDR-2 is 
calibrated by a time signal receiver (tuned to the coded time signal emitted by the DCF-77 
Station in Frankfurt, FRG). The NARS array comprises 14 of these mobile, 3-component, 
digital seismic stations. The stations record in a triggered mode and write their data on 
high density cassette tapes. 

For the NARS experiment in western Europe (Nolet & Vlaar, 1982; Dost et al., 1984; 
Nolet et al., 1985) the seismometers were adjusted to an eigenperiod of 12 s, and the 
shaping filter modified this response to that of a seismometer with a 100 s eigenperiod. The 
sampling period was set to 0.125 s, and a second-order low-pass Butterworth filter with a 
comer frequency of 1 Hz (in the event recorder) was applied as anti-alias filter. The NARS 
instrument response to velocity and the responses of the individual components of this set- 
up are shown in fig. 2.2, The total response to velocity is approximately flat in the 10 to 
1000 mHz frequency band, making it essentially a broad-band’ station. 
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NARS INSTRUMENTATION 


SL 210 
[as CHOPPER PRE- 
ae AMPLIFIER 
SL 220 


PDR-2 
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SEISMOGRAPH 
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TIME FIELD 
SIGNAL |_| DECODER PLAYBACK 


RECEIVER 


Figure 2.1. The NARS system design 


NARS amplitude response to velocity 


Total response 

Anti-alias filter 

Seismometer 
— — — Shaping filter 


(V/ms) 


log gain 


log frequency (Hz) 


Figure 2.2. Instrument response of the NARS stations. The amplitude in V/ms is converted to 
counts/ms by the factor 104857.6 counts/Volt, 
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With the Fourier sign convention 
F(a) = fe dw 
the transfer function to displacement is given by (Dost et al., 1984) 


-iA@ 


T(@) = D(@) 


The denominator D can be expressed as 
D(@) = (—@* — 2imw,h, + 2) (-0? — 2iww,h, + 2) 

with 

@, =2n/100 (rad/s) 

@, =2n (rad/s) 
=1.0 
h, = 0.707 
A =1.23x10'° (counts/ms ) 


a 


Alternatively, the transfer function can be described in the form of a 4-pole recursive filter 
where 


D(@) = (+ W; )( @ + Wy )( @ + 03 (CO + My) 


with poles 
@, = 4443+4.443i (rad/s) 
@ =-4.443+4,443i (rad/s) 
@, = 0.0628: (rad/s) 
@, = 0.0628: (rad/s) 


2.2 STATION CONFIGURATION AND OPERATION OF THE NETWORK 


For the NARS project in western Europe, the stations have been installed in a 
configuration optimal for (higher-mode) surface wave dispersion measurements from 
Japanese earthquakes. The stations were therefore located in 6 different countries 
approximately along a great-circle that passes through the region of high seismic activity 
near Japan and the station spacing was about 200 km. Although the stations were 
preferably installed in observatories, the majority of the sites were installed in less well 
equiped places. In some cases these locations proved to be inadequate, and stations had to 
be replaced. Figure 2.3 shows the configuration of the NARS array, and table 2.1 presents 
the coordinates and the period of operation of the stations. 
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Figure 2.3. NARS station configuration. 


Station 
Gothenborg 
Monsted 
Logumkloster 
Witteveen 
Utrecht 
Dourbes 
Villiers-Adam 
Aigurande 
Les Eyzies 
Arette 
La Almunia 
Ainzon 
Valle de los Caidos 
Puertollano 
Granada 
Valkenburg 
Clermont Ferrand 
Toledo 
Rejaudoux 


57.801N 12.132E 
56.459N 9.170E 
55.045N 9.153E 
52.813N 6.668E 
52.088N 5.172E 
50.097N 4.595E 
49.074N 2.232E 
46.420N 1.730E 
44.852N 0.981E 
43.086N 0.699W 
41.477N 1.372W 
41.814N 1.517W 
40.642N 4.155W 
38.685N 4.091 W 
37.190N 3.595 W 
50.867N 5.785E 
45.763N 3.103E 
39.881N 4.049W 
45.304N 1.516E 


1 Station closed between Jun. 1984 and Jan. 1986 


Period of operation 
Feb. 83 - Jan. 86 
Feb. 83 - Mar. 87 
Feb. 83 - Mar. 87 
Jul. 82 - Feb. 88 
Mar. 82 - Feb. 88! 
Jul. 82 - Feb. 87 
Nov. 83 - Jun. 87 
Nov. 82 - Dec. 84 
Nov. 82 - Feb. 86 
Nov. 82 - May 87 
May 83 - Nov. 83 
Nov. 83 - Dec. 87 
May 83 - Feb. 85 
May 83 

May 83 

Jun. 84 

Nov. 84 - Jun. 87 
Feb. 85 

Feb. 86 - Jun. 87 
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The task of the local station managers (mostly volunteers) primarily consisted of 


Event 1983-08-06, broad-band and short-period seismograms 


2800 


BER-SP 


fl 


L 


NE02-BB 


2600 


2400 


2200 


Distance (km) 


2000 


Time - Distance / 10.0 (8) 


Figure 2.4. Broad-band NARS, and short-period DWWSSN P-wave arrivals on the vertical 
component of an event in the Aegean Sea (Date: 1983-08-06, T,: 15:43:51.9, Coord: 40.14N, 
24.74E, Depth: 2 km, m,: 6.0). 


changing cassette tapes at least once every two weeks and checking the operation of the 
Station. The cassette tapes were sent to Utrecht, where the data were processed and 
converted to a format consistent with that of the GDSN (Global Digital Seismograph 
Network) event tapes. 


2.3 ADEQUACY OF BROAD-BAND SEISMOGRAMS FOR BODY WAVE 
STUDIES 


This thesis presents the results of two types body wave studies of the upper mantle. I 
found that the broad frequency response of the instruments was ideally suited for these 
purposes, enabling P and S phases to be modelled using the same data set (chapter 4), and 
yielding high-frequency information from the upper mantle discontinuities (chapter 5) that 
would not be available from conventional long-period instruments. 

For the study of the upper mantle seismic velocity structure beneath Europe 
seismograms have been used with epicentral distances in the range of 15° to 25°. In this 
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Event 1986-06-18, Station NEO6, broad-band Event 1986-06-18, Station GRFO, long-period 


Delta: 78.6, Azm: 1.0, Depth: 61 Delta: 78.8, Azm: 5.2, Depth: 61 


720 740 760 780 800 720 760 800 840 880 
Time since origin (s) Time since origin (s) 


Figure 2.5. Broad-band (a) and long-period (b) seismograms of an event from the Aleutian 
Islands (Date: 1986-06-18, T,: 08:05:16.0, Coord.: 51.66N 176.96W, Depth: 61 km, m,;: 5.7). 
The three components are marked by Z (vertical), R (radial) and T (transverse). The epicentral 
distance and backward azimuth to the stations NE06 (a), and GRFO (b) are given in the figures. 
Indicated with the arrow in (a) is a possible P-to-S converted phase from the 670-km 
discontinuity; the corresponding location has been indicated by the question mark in figure (b). 


distance range, multiple arrivals from the low velocity zone and sharp velocity gradients 
beneath it produce a complex signal. Identification of the different arrivals is necessary to 
obtain at least the travel time information of the phases. In chapter 4 (section 2, see fig. 3) 
it is shown that an unambiguous interpretation of the S-wave arrivals is not possible on the 
basis of long-period seismograms, and Ingate et al. (1986) demonstrated that this hampers a 
detailed interpretation in terms of the upper mantle velocity structure. On the other hand, 
one might try to model the (higher-frequency) P phases using short-period seismograms 
(maximum amplification at | Hz). Figure 2.4 shows the P-wave arrivals on the vertical 
component of two short-period and three broad-band (NARS) stations. The data are for the 
same event as for the (S-wave) example shown in chapter 4 (fig. 3; P phases have not been 
recorded by station NE12, which is located near to station TOL, because triggering 
occurred on the S-wave arrivals). It is clear from fig. 2.4 that the long-period energy is 
absent on the short-period seismograms. Furthermore, these data show more ’ringing’ with 
a dominant period of about 1.5 s. This is probably due to shallow crustal reverberations 
and wavetype conversions at the near-receiver structure. This high-frequency signal tends 
to mask the upper mantle arrivals; low-pass filtering might in this case facilitate the 
interpretation by enhancement the low-frequency content. 
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Fig 2.5 has been included to show the difference between the broad-band and long- 
period seismograms for the interpretation of P-to-S converted phases from the upper 
mantle. Note the difference in time scale between fig. (a) and (b) by a factor of 2, The 
arrow in fig. 2.5(a) points to the onset of a possible P-to-S conversion from the 670-km 
discontinuity on a seismogram of NARS station NE06. Figure (b) shows the long-period 
seismogram of the same event recorded in station GRFO (with only slightly different 
epicentral distance and backward azimuth compared to NE06). The question mark 
indicates the corresponding location for the P-to-S converted phase on the long-period 
seismogram of GRFO. It is evident that these low-amplitude phases are difficult to identify 
on long-period seismograms, due to the lack of characteristic high-frequency signal from 
the source time function. Moreover, when identified, it does not allow a detailed 
interpretation in terms of the width of the 670-km discontinuity. 

In conclusion, many of the inferences that are presented in these thesis, could not have 
been made if the instrument response of the NARS stations were peaked at a certain 
dominant frequency. 
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Chapter 3 


The WKBJ method 


In the last two decades synthetic seismograms have provided a valuable tool for the 
interpretation and inversion of seismic data. Many different techniques have been 
developed to calculate theoretical seismograms in horizontally stratified media, and, more 
recently, also for laterally inhomogeneous media. For practical reasons, all methods 
involve approximations. The WKBJ method (Chapman, 1978) is a technique where the 
high-frequency seismic signal is approximated for a medium consisting of vertically 
inhomogeneous layers. I have used the method extensively in this thesis, especially in 
chapter 4, because it efficiently and adequately calculates body wave seismograms for a 
large number of signals. In some cases, however, the technique introduces errors that are 
inherent to the approximations involved. This chapter intends to investigate the conditions 
for its validity and the cases where it breaks down. For a more elaborate description of the 
WKB] seismogram method, particularly for the derivation of the elastodynamic solution, 
one should refer to Chapman (1978), Dey-Sarkar & Chapman (1978), or the review article 
by Chapman & Orcutt (1985). 


3.1 WKBJ THEORY 


In its most general sense, the WKB theory provides an approximate solution to a 
linear differential equation whose highest derivative is multiplied by a small parameter. It 
is named after Wentzel, Kramer, and Brillouin, who rediscovered the theory and made it 
accessible for general use. In seismology, it is generally called WKBJ theory, because 
Jeffreys (but also Rayleigh) contributed to its early development. The WKBJ theory in 
seismology is concerned with the solution of the wave equation for high frequencies, and in 
the following, an outline is presented for the solution of a general (2-D) scalar wave 
equation in the frequency-space domain. Its generalization to the frequency-slowness 
domain (the WKB] seismogram) is briefly summarized in the next section. 
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We start from the 2-D (x = (x, 2)") scalar wave equation in the frequency domain: 
V7o + w? v-(2) 6 =0 (1) 


The physical meaning of o and v(z) is not important at this point, but these variables could 
be interpreted as the pressure and the acoustic velocity, respectively. To find the 
approximate solutions for equation (1) from WKBJ theory, it is required that o is large. 
Trial solutions are of the form 


ad i 
(0, x)= 9,(0) EAM gia ©) 
n=o (-I0) 
Substitution of this expansion in (1) and requiring that the coefficients of terms of order w* 
sum to zero yields the eikonal equation 


(VT) - v-(z) = 0 (3) 


We define p = VT for real solutions of VT. The slowness vector, p = (p, 4)‘, is 
perpendicular to ’wavefronts’ of constant ¢ = T(x). The amplitude coefficients A(x) can 
be obtained by setting the sum of terms of each power of w equal to zero. 

At high frequencies, only the leading term of the WKBJ asymptotic expansion will be 
important, giving the WKBJ approximation 


(2, x) = o,(@) AO (x) e7™™ (4) 
or in the time domain 
(t, x) = A(x) o,[¢ — T(x)] (5) 


the well-known solution from geometrical ray theory for the propagation of a pulse o, with 
amplitude variation (’ geometrical spreading’) A(x) and travel time T(x). 

For the two-point boundary value problem, it will be found that often several solutions 
for T(x) and A™(x) exist. Therefore, in the zeroth-order WKBJ approximation, we can 
write 


ra. 
$(0, x) = 4,(0) An) ) 
fl 
where ’rays’ represents the number of solutions. 

It is important to note that the zeroth-order approximation is only valid if the velocity 
distribution is continuous. Or, more generally, in order to determine the nth term in the 
asymptotic expansion, it is necessary that the nth derivative of the model parameter is 
continuous. Thus, the WKBJ approximation breaks down if a discontinuity in v(z) is 
encountered. The solution is then propagated across the interface using Snell’s law and the 
appropriate reflection and transmission coefficients. 

Furthermore, it is essential to understand that we actually have two families of non- 
decaying solutions, as from equation (3) 
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p+q =v) 
or 


q(p, 2) = £Vv(z) ~ p? (7) 


Waves traveling downward (g > 0) and waves traveling upward (g < 0) are coupled at g(p, 
z) = 0, where p = v!(z). However, exactly at the branch point for g = 0, the turning point 
of the ray, the WKBJ expansion breaks down. This problem can mathematically be 
overcome by connecting the oscillatory behaviour above the turning point (real phase 
function) to the exponentially decaying behaviour below it (imaginary phase function) by 
locally applying the Langer approximation (see e.g. Chapman, 1974), 


3.2 THE WKBJ SEISMOGRAM 
GENERALIZATION OF GEOMETRICAL RAY THEORY 


The generalization of the geometrical ray theory described in the previous section to 
the WKBJ seismogram involves an extra Fourier transform of one of the spatial coordinates 
to the slowness domain. The essential point of the WKBJ seismogram is the order of the 
back transformation: the solution in the frequency-slowness domain is first transformed to 
the time domain before the transformation to the spatial coordinate is applied. 

The Fourier transform with respect to the x coordinate is defined as: 


We 
6(@, p, z) = -2 J (a, x, z) eo * dx (8a) 
and the inverse: 
20 
(0, x, z) = | io] J (0, p, z) &* dp (8b) 


The zeroth-order term of the WKBJ asymptotic expansion in the frequency-slowness 
domain can be written as 


rays 5 
(0, P, 2) = bo) 1 AS%(x) eH” (9) 
Fl 
Using the transformation (8a), it can be shown that the phase functions T(x) of eq. (6) and 
t(p, z) of eq. (9) for a certain slowness p are related by (Dey-Sarkar & Chapman, 1978) 
Up, z) = T(x(p), 2) — x(p) p (10) 


The inverse transform to the time domain can be solved analytically, and the inverse 
slowness transform then becomes (Dey-Sarkar & Chapman, 1978) 
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(4, ») =e 4,(0 * 3, Im} AG) * [ AG, 2) 8(¢- 0, x) ) dp (11) 
m2 “ve 


where 
AQD=MDN +i MO 


with A(t)=H(t) "2 and 4) =H(-1) (-)"?, the Hilbert transform of X(t), H(t) the 
Heaviside function, and 


O(p, x) = t(p, z) + px 
Contributions to the integrand of (11) are obtained for ¢ = 8(p, x), giving 


a (0) 

WO) =e 0,08, in| 0 : Be (12) 
For all real ¢ = 6, the inverse transform is exact. However, for some slownesses p, T(p, z) is 
complex, producing an exponentially decaying phase term (with a high decay for high «). 
In the WKBJ seismogram method the inverse transform is evaluated only over a real p 
contour. In practice, the (real) slowness range is limited either by numerical truncation or 
by complex slowness values. This discontinuous behaviour of the amplitude terms causes 
end-point errors on the WKBJ seismograms, as will be shown in the next paragraph. 

It should be noted that the most important contributions to the solution of eq. (12) are 
obtained for ¢ = ® with 0, 0, while singularities exist for 0, = 0. The physical 
interpretation of these results as applied to the seismological problem will be discussed in 
the following section. 


SOLUTION FOR AN ELASTIC, LAYERED MEDIUM WITH A SOURCE TERM 


Expressions (11) and (12) describe the wave propagation through an inhomogeneous 
medium without discontinuities, the solution of (1). To obtain the solution for more 
realistic seismological problems, we must calculate the displacement due to a point force in 
an (isotropic) elastic medium in which discontinuities are present. Solutions can be 
obtained from the elastic potentials ( in eq. (1) can interpreted as the P-wave potential: 
o = V-u, , where u, represents the displacement), but when the models contains interfaces it 
is more convenient to represent and solve for the components of displacement and stress in 
the form of a vector, the displacement-stress vector. The zeroth-order term of the WKBJ 
expansion is then obtained in the same way as described in the foregoing, but the solution 
in terms of the displacement-stress vector involves more ’bookkeeping’. For laterally 
homogeneous media, we find the well-known results of P waves that propagate with a 
velocity o(z) and that have a longitudinal polarization, and S waves, that can be separated 
into two perpendicular transverse SV and SH components, propagating with a velocity B(z). 
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Interfaces are now easily incorporated by including the appropriate reflection and 
transmission coefficients of the rays. By the requirement of continuity of stress and 
displacement at each interface, the P and SV components combine to form a coupled 
system. The introduction of a point force f =f,(t) 6(x—x,), does not present problems in 
the WKBJ approximation. The solution of the now inhomogeneous wave equation (with a 
term V-p'f, on the r.h.s. of eq. (1)) is easily obtained at x,, and the general solution is 
obtained by propagating the solution to the receiver at x. A moment tensor instead of a 
point force as source term is represented as m= V- M,,(f) 8(x—x,) (Chapman & Orcutt, 
1985). 

Following this procedure, and expanding the displacements and stresses to vector 
cylindrical harmonics, Chapman (1978) and Dey-Sarkar & Chapman (1978) showed that 
(to zeroth order) the WKBJ seismogram for P waves at x = (x, z)' due to an explosive point 
force is 


n(t,x)=- — FO _ em) aw * f LR 4 8-8) dp] (13) 
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where u represents the displacement, fe) the time derivative of the point force f(f) at x, = 
(0, 0)", and the subscript o refers to the model parameters at the source. R(p) is the product 
of all reflection and transmission coefficients along the ray and includes factors of -i sgn(w) 
(zeroth-order Langer approximation) for each turning point. The slowness integral is 
evaluated at 


O(p, x)= t(p) + px (14) 


with phase function 
1(p) = J ap, 2’) dz’ (15) 


where the integral from z, to z represents the ray integral over z between source and 
receiver. The travel time and phase function are related by (analogous to eq. (11)) 


(p) = T(p) — x(p) p (16) 
The most important contribution to the integral in equation (13), for a slowly varying 
function R(p), is obtained for 0,0(p,x)—>0 (see fig. 3.1(d-e)), or equivalently, for 
v(p)+ x — 0 (see fig. 3.1(b)). So, if we define X(p) = — 1’(p), then we obtain saddle points 
for p = p, with x = X(p,). It is easily shown that the travel time function T(p,) is equal to 
the theta function at the saddle point p, (as indicated in fig. 3.1(d-e)). To remove the 
singularities in (13) for each ray geometrical arrival where 0,0=0, ¢ is expanded around the 
saddle point p,: 
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Figure 3.1. Travel time curve (a), T-p curve (b), x-p curve (c), 8-p curves for two different 
distances (d and e), and raypaths (f) for an upper mantle model with a triplication. The letters A- 
D indicate different branches or points on the travel time curve (a), and the corresponding points 
have been indicated on the other diagrams. Fig. (d) shows the 6-p curve for the distance at C. 
The travel times (at 0,0 = 0) are indicated by T, and T>, and the time corresponding to the 
maximum slowness by t,. Fig. (e) shows the 0-p curve for a distance of 1805 km. The three ray 


geometrical travel times are indicated by T,—T3 (point C does not coincide with the curve at T4). 


t= O(p, x) = T(p,) + 1/2 (p— p,)* 020 
= T(p,) — 1/2 (p~ p,)’ Xp) (17a) 


For t > T(p,) expression (13) can be evaluated for the two (real) roots p; (i = 1, 2) on either 
side of the saddle point, and the argument of the delta function can be approximated by 


t-O@= (p — pi) 0,9(0;, x) 
= (p— pj) [x-X@) ] (17b) 


Both roots contribute half of the geometrical amplitude of the solution at ¢ = 0. 
EXAMPLE TO SHOW THE RELATIONS BETWEEN p, T, X, 8, AND T 


Figure 3.1 illustrates the relations between the different parameters for a triplication in 
the travel time curve as shown in (a). The t-p curve (b) is directly calculated from integral 
(15), and classical ray theory provides the well-known relationship between x and p ina 
spherical or flat layered earth (fig. c). The distance X for the ray geometrical arrival with a 
slowness p can be determined from this diagram as is indicated. Figures (d) and (e) present 
0-p curves (eq. (14)) for two different values of X. The travel times T,(p,) for different 
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Figure 3.1. (continued) 


branches of the travel time curve are indicated in figures (d) and (e). For completeness the 
raypaths through the model are shown in fig. (f). 
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SMOOTHING 


A straightforward way to approximate the singularities in expression (13) for the 
numerical computation of the seismogram is to smooth the signal in the time domain. The 
program used in chapter 4 convolves the time series with a boxcar window, B(t), of length 
2At, where At is the sample interval. 


Bi) = 1/2 | Ht+1)-AHG¢-1)] 


where H(t) is the Heaviside function. The smoothed solution becomes 


B(t/At) 0) 1 p'? R(p) 4 
xy* 2 ee — * Im] A() * —— ee dp| (17 
Bee ae QP gy? OO” oat poets (ppg)? \e) ” ms 


Thus, the integral is numerically evaluated over finite (real) slowness intervals 
corresponding to t= 6(p,x) + At. Smoothing has the important advantage that the results 
are not sensitive to small features in the model, allowing the model parameters to be 
linearly interpolated without generating unwanted effects. 


EARTH FLATTENING TRANSFORMATION 


Although the WKBJ seismogram approach can be calculated for a medium with 
spherical symmetry as well as for a flat layered medium, it is convenient to use the same 
algorithm for both problems. The transformation of the spherical system to a plane layered 
system involves a simple transformation of the model parameters (see e.g. Aki & Richards, 
1980, p. 463-465). It has been shown that this ’earth flattening transformation’ is exact for 
the SH problem, but involves first order errors (O(w~')) for the P-SV problem; these are of 
the same order as for the WKBJ approximation. 


3.3 EXAMPLES 


In the previous sections it has been pointed out that the WKBJ seismogram is a 
zeroth-order, high-frequency approximation, for which the inverse transform is evaluated 
over a (limited) real p contour, and that the signal is smoothed in the time domain. The 
earth flattening transformation does not introduce extra errors to this approximate solution. 
In this section, the adequacy of the WKBJ method and its shortcomings will be illustrated 
by synthetic seismograms. 


TURNING RAYS CALCULATED FOR A LIMITED SLOWNESS INTERVAL 


The WKBJ method accurately (and economically) calculates synthetic seismograms 
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Figure 3.2, WKBJ seismograms of direct P phases with turning points in the lower mantle of 
model 1066B. End-point arrivals due to truncation of the slowness interval are indicated by the 
arrows, in fig. (a) for the maximum slowness and in fig. (b) for the minimum slowness. 


for ’turning’ rays. Expression (17) adequately describes the signal at t = T(p,), the time 
corresponding to the ray geometrical arrival at a distance x, for a smoothly varying function 
R(p). However, so-called ’end-point arrivals’ are also generated due to the evaluation of 
the slowness integral over a (numerically defined) limited slowness range. The 
discontinuous behaviour of the integrand in (17) at the end-point slownesses p, causes 
end-point arrivals at 


t, = O(p,, x) =T(p.) + p[x-X(p.) ] 


It can be shown that the amplitude of these arrivals decays as 1 /!x — X(p,)| from X(p,) for 
t, > T(p,). Figure 3.2 presents the WKBJ seismograms for turning P waves in the lower 
mantle of model 1066B (Gilbert & Dziewonski, 1975). The time axes of fig. 3.2 (a) and (b) 
have been reduced with reduction velocities corresponding to the maximum and minimum 
slowness of the integrand, which makes the end-point arrivals line up vertically on the 
diagrams. These end-point errors can be confusing, but, in general, their amplitude is 
significantly smaller than that of the ray geometrical arrivals. Moreover, if the difference 
between the ray geometrical slowness and the end-point slowness is sufficiently large, then 
the end-point phase will arrive outside the time interval of interest. However, there may be 
cases that its effect is not negligible. A detailed interpretation of the synthetic seismograms 
can then be misleading if one is not aware of these (predictable) artefacts. 


TRIPLICATIONS OF THE TRAVEL TIME CURVE 


Model 1066B has been used to compute synthetic seismograms of upper mantle P 
phases shown in fig. 3.3. The seismograms drawn with solid lines represent the synthetics 
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Figure 3.3. Synthetic seismograms of upper mantle P phases for model 1066B. The solid line 
marked with A-F indicates the travel time curve. The other solid lines represent the seismograms 
for the original model with first order discontinuities. The dashed lines represent the 
seismograms for the same model where the discontinuities have been replaced by velocity 
gradients extending over a width of 1 km. 


of the original model with discontinuities in the velocity distribution at 420 and 671 km 
depth (’first order discontinuities’), The dashed lines represent the synthetics for the same 
model except that the discontinuities have been replaced by velocity gradients extending 
over a width of 1 km. Thus, for the seismograms computed for the original model it has 
been necessary to propagate the solution across the interfaces using the appropriate 
reflection and transmission coefficients, whereas for the ’gradient model’ the seismograms 
could be calculated using the WKBJ approximation only. It is clear that both sets of 
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seismograms are very similar except for a few features, and these features point to essential 
aspects of the approximations in the WKBJ seismogram. 

It should be noted that, although the WKBJ seismogram is a high-frequency 
approximation, it does not predict infinite amplitudes at caustics. This is a clear 
improvement over geometrical ray theory that gives singularities in the amplitudes at the 
points C and E of the gradient model due to a locally infinite geometrical spreading 
((dxidp)” = 0). However, the WKBJ seismogram does not provide a completely adequate 
solution for velocity gradients, because it does not model low-frequency reflections from 
gradient zones. This is inherent to the WKBJ approximation because the frequency and 
space dependence are separated in the solution. The relative amplitude of the solution at a 
depth z depends only on the model parameters at that depth (see, for instance, eq. (9)). 
Low-frequency reflections from gradient zones are therefore not modelled, even if higher- 
order terms are included in the solution. A way to overcome this problem practically is to 
model the gradient zone by a stack of thin layers and calculate the reflections from each 
interface (and include internally reflected phases for a better approximation). However, 
this is not an efficient method and preference should be given to an alternative technique 
(e.g. *full wave method’, Richards, 1973) to model reflected waves from gradient zones 
adequately at low frequencies. The difference between ’gradient reflections’ and (ray 
theoretical) reflections from a first order discontinuity can be seen in the regions of pre- 
critical reflections (e.g. at x < 1800 km for the 671-km transition). The amplitude of the 
pre-critical ‘gradient reflection’ decays faster than that of the first order discontinuity 
because of its (larger) geometrical spreading. 

Another difference between the two sets of synthetic seismograms is that the 
diffracted phases from the 420-km discontinuity (the extension from B in the 2600 to 3500 
km distance range) have a significantly larger amplitude for the sharp discontinuity than for 
the 1 km thick gradient. This observation points to one of the most important shortcomings 
of the WKBJ seismogram method: the interaction of waves with interfaces at finite 
frequencies is not properly accounted for. The solution for arrivals in the shadow zone is 
obtained by expanding the function 6 (eq. (17a)) around the saddle point p,, the slowness 
corresponding to the grazing ray. The amplitude is then half of the geometrical value at the 
shadow edge, because one of the roots p; is imaginary. For arrivals in the (deep) shadow, 
the WKBJ method just models the effects of Fresnel diffraction (using eq. (17b), yielding 
an amplitude decay of 1 /lx— X(p,)!). The method does not model the interaction of the 
wavefront with the interface at finite frequencies. The amplitudes of diffracted phases in 
the shadow zone for the model with the sharp discontinuities are therefore too large. 
Similar effects are expected for waves turning just above an interface, where the zeroth- 
order Langer approximation (a factor -i sgn(w)) is used to model the signal, but wave 
interaction with the discontinuity is not negligible. 

The slightly larger amplitudes of the arrivals for the gradient model on the branch EF 
(turning phases below the 671-km transition) are most likely due to the fact that more 
energy is transmitted downward from shallower depths for this model than for the model 
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with the interfaces. These differences are not due to approximations in the method. 

It can be concluded that, when modelling arrivals from a velocity structure with 
interfaces, one should be aware that amplitudes of diffracted phases are overestimated. A 
for most purposes less inconvenient shortcoming is that low-frequency energy reflected 
from gradient zones is not modelled by the WKBJ method. 


LID AND LOW VELOCITY ZONE 


Lid and low velocity zone 
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Figure 3.4. (a) Solid lines represent the WKBJ seismograms of P phases of the velocity structure 
shown in (b). The short dashed lines indicate the travel time curve. The thick long dashed lines 
mark the diffracted phases, and the thin long dashed lines mark end point arrivals. Fig. (c) shows 
the raypaths through the model. 


Figure 3.4(a) shows the WKBJ seismograms of direct P phases for a model with a low 
velocity zone. The model and raypaths through the medium are presented in fig. 3.4 (b) 
and (c), respectively. Ray geometrical arrivals from the high velocity ’lid’ (branch A) and 
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Figure 3.4. (Continued) 


from the velocity structure below it (branch BC and CD) are easily recognized on fig. 
3.4(a). Diffracted phases in the shadow zone from the caustic at B and from the lid (A) are 
also clear, whereas end-point arrivals evidently contaminate the signal after the first arrival 
in the 600-1800 km distance range. It has already been pointed out that the WKBJ method 
adequately calculates the high-frequency signal at caustics, but only models the effect of 
Fresnel diffraction in the shadow zone. These aspects have been inferred for the case of a 
triplication, but the same arguments apply for the shadow of a low velocity zone (at C). 
Although it cannot be inferred from the seismograms displayed in fig. 3.4(a), it is 
expected that the amplitude of the diffracted phase from the base of the lid is too small, 
since the WKBJ method does not model the tunneling of low frequency energy at the base 
of the lid. This is, in essence, the same shortcoming as that for gradient zones discussed in 
the previous example. However, in this case it is not possible to approximate the low- 
frequency effects by adjusting the model. The WKBJ method is therefore inappropriate for 
calculating (non-ray-geometrical) signals from a lid-low-velocity structure. These phases 
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may, however, be important, especially for modelling upper mantle phases for a velocity 
distribution of a tectonic region (see e.g. chapter 4, fig. 15(a): the Pn phase at NEO2). One 
should in such cases be aware of the shortcomings of the WKBJ technique, and not try to 
model these low-frequency arrivals, or use a more adequate technique for these problems. 


3.4 CONCLUSIONS 


Several examples of WKBJ seismograms for structures relevant to the upper mantle 
have been discussed in the previous section. It was shown that ray geometrical arrivals are 
properly calculated, but low-frequency effects from velocity gradients and wave interaction 
with interfaces are not modelled. The adequacy of the method for the majority of the 
arrivals, combined with the effectiveness of the computations (although not mentioned 
explicitly) makes this method a suitable technique for ’reconnaissance studies’ of 
waveform data that are inverted to a velocity structure. These aspects, together with the 
fact that its shortcomings and artefacts are easily predicted, have led to the choice of this 
technique for the study of the upper mantle beneath Europe as presented in the next 
chapter. It should be emphasized that the technique has been used in the next chapter 
mainly as an aid to interpret the multiple upper mantle arrivals and to match their travel 
times. More detailed (amplitude) modelling is not justified considering the shortcomings of 
the WKBJ method, particularly for diffracted phases from low velocity zones and 
discontinuities. 
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Chapter 4 


Lateral heterogeneity of Europe’s upper mantle 


Summary. The upper mantle velocity structure beneath Europe is investigated 
by body-wave modelling of broad-band seismograms from earthquakes in the 
eastern Mediterranean recorded by stations in north-west Europe. The gross 
features of the P- and S-wave velocity distribution as inferred from the data 
show a high velocity lid, a low velocity zone, and two discontinuities at 
depths of approximately 420 and 650km. The seismograms also present 
evidence for lateral heterogeneity of at least a few per cent (2.5 per cent for 
P, 4 per cent for S) to a depth of approximately 400 km. Many of the 
inferred velocity variations can be related to the tectonic setting of the region 
between source and receiver. Broad-band seismograms are essential for such 
detailed study of upper mantle structure using body-wave modelling. 


1 Introduction 


In recent years many studies have been carried out to investigate the lateral heterogeneity of 
the Earth’s upper mantle. The objective is to gain better insight into the geodynamical 
processes occurring in the mantle and to obtain a better knowledge of its chemical 
composition. At present, large-scale heterogeneity in the seismic velocities up to a few per 
cent has been established to depths of about 400 km (e.g. Woodhouse & Dziewonski 1984; 
Grand & Helmberger 1985). But, in spite of the results produced by various studies, it is still 
impossible to discriminate between different concepts of the upper mantle such as suggested 
by Ringwood (1975), Anderson (1979, 1984), Jordan (1978), or Vlaar (1983) (see Nolet, 
Dost & Paulssen 1986 for a more detailed discussion on this subject). One of the important 
questions in this respect is to what (small) scale distinct lateral heterogeneity may be 
present. 

Methods which have been used to investigate the upper mantle structure include funda- 
mental and higher mode surface-wave analysis. These data have yielded the gross S-wave 
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velocity structure to a few hundred kilometres depth (e.g. Woodhouse & Dziewonski 1984; 
Nolet 1977; Cara, Nercessian & Nolet 1980; Dost 1986), but lack the resolution on a 
smaller scale, say less than 1000 km horizontally and 50 km vertically. Tomography, using 
P arrival times, points to lateral variations in the P-wave velocity structure in Europe which 
could be established within a horizontal length scale of 100 km (Spakman 1986a). However, 
the gradients of the velocity—depth function, and the depths of the discontinuities cannot 
be determined accurately with tomographic studies. This information is essential for a 
petrological interpretation of the upper mantle velocity variations. Relatively accurate 
velocity models are abtained by body-wave analysis of regional data, and many models 
representing the average upper mantle velocity structure of specific regions have been 
determined in the last decade (e.g. Burdick 1981; Walck 1984, 1985; Leven 1985). The 
variations among the models suggest a large degree of lateral heterogeneity, but the 
horizontal length scale on which this occurs cannot be estimated. 

This paper presents the results of a waveform-fil analysis of upper mantle phases with ray- 
paths in Europe. [t is not the first study which adopts the technique of body-wave modelling 
for a determination of the upper-mantle velocity structure, but there is an important 
difference with previous studies of this kind: the use of broad-band seismograms. Most 
previous studies have employed long-period data to infer the average velocity structure of 
the area (e.g. Helmberger & Engen 1974; Burdick & Helmberger 1978: Lyon-Caen 1986). As 
Ingate, Ha & Muirhead (1986) have pointed out, there is a severe limitation on the resolution 
which can be obtained by working with long-period seismograms. Only the gross features of 
the velocity structure can be resolved. The results of this study show that arrivals from 
strong velocity gradients and low velocity zones can be identified more accurately on broad- 
band seismograms, thus enhancing the resolution of the inferred velocity model. 

The data which have been analysed are from earthquakes in the eastern Mediterranean 
recorded by NARS stations in north-west Europe (NARS: Network of Autonomously 
Recording Seismographs, cf. Nolet & Vlaar 1982, Dost, van Wettum & Nolet 1984). P- and 
SH-wavetorms are modelled by trial and error using the WKBJ technique (Chapman 1978) to 
yield the upper mantle P- and S-wave velocity structure. The waveform fitting results show 
that the average velocity structure for the region between south-east and north-west Europe 
has a high velocity lid, a low velocity zone, and two discontinuities at depths of approxi- 
mately 420 and 650 km. The data, however, also present evidence for lateral heterogeneity: 
it is not possible to model the seismograms with a single P- or S-wave velocity distribution, 
and different models must be assumed to match the arrival times of the phases on different 
seismograms. The degree of lateral heterogeneity of the area is roughly estimated by 
comparison of 1-D models which are thought representative of the average structure along 
the ray-paths. The results suggest that velocity variations of at least a few per cent to the 
depth of the 400-km discontinuity are present in the upper mantle beneath central Europe. 

Al present, it is not possible to perform a 3-D inversion by waveform analysis because of 
the sparseness of the available data. However, the results of this study indicate that 
modelling of broad-band seismograms, which are now becoming more widely available, 
allows a more detailed determination of upper mantle structure. 


2 Data 


The data which have been used in this study are both P- and SH-waves for events at distances 
of 15°-25°, recorded by broad-band stations of the NARS array. Events and stations are 
shown in Fig. 1, and Table | presents the event parameters as published by ISC, or PDE 
when the ISC data were not yet available (event 5-270). All events in the eastern 
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Figure 1. Locations of NARS stations and eastern Mediterranean events. 


Mediterranean which triggered the recorders on body wave phases have been used. The study 
was restricted to the seismograms of stations NEO1—NE16, because these data provided the 
densest sampling of the upper mantle between source and receiver region, thus allowing a 
detailed investigation of this part of the upper mantle. 

An example of P- and SH-wave data for an event near Crete (event 4-173) is shown in 
Fig. 2. The observed variations in waveform of the wavetrains can be attributed to several 
effects. First, changes in epicentral distance can easily lead to large variations in the wave- 
forms. Especially in the range of interest (15°—25°) interference of reflected and refracted 
waves can be severe. The radiation pattern of the source can also cause large changes in wave- 


Table 1. Event parameters. 
code date _ origin time coordinates depth m, 
2-229 «17 81982 = 22:22:20.0 33.71N 22.94E 1 6.1 


3-017. 17: 11983 12:41:30.1 38.07N 20.25E 14 6.2 
3-082. 23: 3.1983 = 23:51:05.5 38.23N 20.29E 13. 5.6 


3-186 § 71983 12:01:27.0 40.33N 27.12E 7 5.5 
3-218 6 81983 = 15:43:51.9 40.14N 24.74E 2 6.0 
4-042 11 21984 8:02:50.5  38.38N 22.08E 19 5.3 


4-173, 21 61984 =: 10:43:40.5 35.31N 23.28E 25 (5.8 
5-270 2791985 16:39:48.7 34.51N 26.60E 61 5.6 
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Figure 2, P-wave data on the vertical component (a), and S-wave data on the transverse component (b) of 
event 4-173, The reduced travel time in seconds is marked along the horizontal axis. 


form signature, particularly for large differences in source azimuth. Finally, the structure 
below the stations may have a large effect on the recorded signals. It is of great importance 
that all these effects are correctly modelled or taken into account in order to obtain an 
unbiased estimate of the lateral heterogeneity of the region under study. A further 
discussion on this subject is given in Section 4. 

Most previous studies of body-wave modelling have used long-period seismograms. As 
Ingate et al. (1986) have shown, only the gross velocity structure can be resolved from this 
kind of data. Fig. 3(a) shows four broad-band NARS and two long-period DWWSSN seismo- 
grams of an event in the Aegean Sea (event 3-218, see Table 1). The difference in frequency 
content of the seismograms of NE12 (NARS) and TOL (DWWSSN) can be seen from Fig. 
3(b). The time separation for different branches of the S-wave travel-time curve varies 
typically between O and 20 s for epicentral distances between 15° and 30°. Thus, it is clear 
that separate arrivals generally cannot be recognized on long-period seismograms (see Fig. 
3a). In contrast, broad-band seismograms, obtained from instruments with a flat frequency 
response to velocity between 0.01 and 1.0 Hz such as NARS instruments, are suitable for 
analysis of the wavetrain of upper mantle phases. Using this kind of data, travel times of 
most upper mantle phases can be determined relatively accurately and easily by synthetic 
modelling of the arrivals. Therefore, it is expected that a higher resolution of the upper 
mantle velocity structure can be obtained. 


3 Method 


In studies of body-wave modelling, the upper mantle velocity structure is determined by 
matching the wavetrains of upper mantle phases. The advantage over studies which use 
travel times of the first arrivals only is that travel-time information of later arrivals is also 
used, as well as amplitude information. In many cases the slowness of first and later arrivals 
can be determined by correlation of seismograms with different epicentral distances. 
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Figure 3. (a) Broad-band NARS, and long-period DWWSSN S-wave data on the transverse component of 
event 3-218. The waveshape of the DWWSSN long-period data is affected by clipping of the instrument. 
(b) Frequency content of the seismograms of NE12 (above) and TOL (below) for the time window shown 


in Fig. 3(a). 


Practically, travel times are the easiest to model, although frequently one will not be able to 
identify the phases without the aid of synthetic seismograms and their travel-time curves. 
The travel-time branches of the direct (not surface reflected) waves are therefore indicated 


on the figures of Sections 5 and 6 which show data and synthetics. 
The data are modelled by trial and error using the WKBJ method (Chapman 1978). 


Since the program we have available calculates seismograms for laterally homogeneous earth 
models it may seem inappropriate to use this technique for this application. However, 
because the horizontal! velocity variations are likely to be several times smaller than the 
vertical variations (the data which are used sample a narrow region), it is assumed that 
effects due to lateral heterogeneity are of secondary importance, and need not be taken 
into account in a ‘first-order approximation’. The method will therefore yield average 
models along the ray-paths, and the degree of heterogeneity is estimated by differences 
between the models obtained for different ray-paths, as far as they cannot be reconciled in a 


single model. 
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Features that are present in the data which cannot be modelled by the WKBJ method 
may sometimes be caused by shortcomings of the WKBJ technique, but the cases where the 
WKBJ method fails are known, and one can deal with these situations. For instance, it is 
known that diffracted rays are modelled with an amplitude which is slightly too large, it is 
awkward to model the low-frequency reflections from strong gradient zones, and tunnelling 
waves cannot be included in the synthetics (Grand & Helmberger 1984: Chapman & Orcutt 
1985). 

Using the WKBJ algorithm, Green’s functions are computed for the P-, pP, and sP-phases 
to model the P-wavetrains on the vertical component, and for S- and sS-phases to model the 
S-wavetrains on the transverse component. The radiation pattern given by a double couple 
solution is incorporated in the synthetics by multiplying the functions of each phase with 
the appropriate weights. For each event a clear P- and S-pulse was selected from the data. It 
is assumed that these wavelets represent the total effect of source time function, instrument 
response and attenuation for the P- and S-phases respectively. Green’s functions are con- 
volved with these pulses to produce the synthetic seismograms. This means that it is 
implicitly assumed that attenuation is the same for all S phases (including sS). and for all 
P phases (including pP and sP). This assumption is not ideal, but upper-mantle attenuation 
is not known accurately enough to be modelled with a high reliability. The results show that 
the approximation is adequate for modelling the most pronounced features of the seismo- 
grams. The consequences for the resolution of the models at discontinuities are discussed at 
the end of Sections 5 and 6. 

Receiver structure effects have not been modelled because they make a relatively small 
contribution to the data. The first Moho reverberation near the receiver has a maximum 
amplitude of about 17 per cent of the direct wave tor the most extreme case: reflections 
within the SH-system and a large velocity contrast across the Moho (3.40--4.76 km 7). 
Signals of this amplitude which arrive in the tail of the wavetrain generally cannot be 
discerned. Constructive interference of Moho reverberations, generating a channelled crustal 
wave, may result in a high-amplitude, long-period wavetrain, but this phase arrives in the tail 
or outside the time window of interest. 

Thus, only the most important phases which may be present in the data are modelled in 
this study. At this stage it does not seem justified to model al! details of the seismograms 
since the structure of the Earth is not known accurately enough to estimate the relative 
importance of the many other contributions. As we will see, the modelling results appear to 
justify this assumption. 


4 The influence of source parameters in studies of body-wave modelling 


It is important to estimate the effects which variations of the source parameters have on the 
synthetically calculated seismograms. The influence of the different parameters on the 
synthetics is illustrated by modelling four of the seismograms (SH-wave data) of the Crete 
event (4-173). 

Synthetic seismograms of the S- and sS-phases on the transverse component calculated for 
the event parameters as published in the ISC bulletin (see Table 1), and in the PDE monthly 
bulletin (To: 10:43:42.1, event coord.: 35.360N, 23.238E: focal depth: 39 km: Fig. 4a) 
show the influence of the event location on the epicentral distance. The focal depth has its 
most pronounced effect on the (relative) travel times of the phases. The difference in origin 
time between the PDE and ISC determinations (1.6 s) is expressed by a shift in the travel 
times of the data which cannot be shown in this figure. Although the differences between 
the ISC and PDE locations for this earthquake were the largest of the events studied, these 
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Figure 4, (a) Synthetic seismograms of event 4-173 calculated for ISC (solid line) and PDE (dashed line) 
event parameters. (b) Synthetic seismograms of event 4-173 calculated for the fault plane solution (solid 
line), and centroid moment tensor (dashed line). 


results clearly show that one has to be careful not to use bulletin locations without further 
consideration. Differences in the event parameters such as these for the Crete event can 
already lead to significant changes in the derived velocity model. Therefore, it is necessary to 
discuss the effects of errors on the modelling results. The differences between the ISC and 
PDE locations and origin times have been taken as an estimate of uncertainty. Travel times 
and synthetics have been calculated for both sets, and it is observed that the arrival times of 
the first onsets differs by 0.5 s or less, which is approximately equal to the accuracy to 
which the P-wave data have been modelled, and much smaller than the differences in origin 
time. This is evidently due to the dependence in the origin time calculations on source 
location, especially focal depth. The event depth affects the time span between the direct 
and surface reflected phases, but for most shallow earthquakes it is not possible to identify 
the pP, sP, and sS phases without the aid of synthetic seismograms. It is observed that the 
synthetics for ISC event parameters given an equal or better fit to the data than the ‘PDE 
synthetics’. A correct estimation of the source depth appears most important for 10-15 km 
deep events, for which early surface reflected phases may interfere with later arriving direct 
phases. The 1SC event depths of the earthquakes used in this study were thought accurate 
enough because of the good fit of the S—sS- and P—pP—sP-waveforms. 

Knowledge of the moment tensor is also of great importance in studies of waveform 
modelling. Fig. 4(b) shows the synthetic seismograms computed for the fault plane solution 
(‘moderately well controlled’), and the best double couple solution of the (Harvard) centroid 
moment tensor as published in the JISC and PDE bulletins. Moment tensor solutions for 
smaller events (m, < 6.0) are generally not well constrained, and although differences 
between the fault plane solution, moment tensor solution and centroid moment tensor may 
be expected, due to the different approaches for the determination, completely different 
(best) double-couple solutions are frequently observed, sometimes resulting in polarity 
reversals of phases over a broad range of azimuths (see Fig. 4b). The moment tensor solution 
has a large effect on the sS/S (or pP/P and sP/P) amplitude ratio on a single seismogram. 
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Thus, good knowledge of the radiation pattern is needed for a good amplitude fit of the 
phases which leave the source at different take-off angles. It is observed that the centroid 
moment tensor gives the best fit to the data for most of the earthquakes of this study, and 
this solution was used to model the data. In view of the uncertainties in the radiation pattern 
and other effects which might influence the amplitudes, modelling efforts were mainly 
concentrated on matching the travel times of the phases. 


5 Determination of the S-wave velocity structure 


The data which have been analysed cover a narrow region as can be seen from Fig. |. The 
best fitting models to the seismograms are found by trial and error. Model 1066B (Gilbert & 
Dziewonski 1975) proved to be an adequate initial (reference) earth model for fitting the 
S-wave data. The travel times of the first arrivals of this model generally correspond within 
2 or 3s with the data. These travei-time anomalies could not be attributed to erroneous 
origin time determinations only, since there appeared to be no systematic positive or 
negative delay for each event. Also, no general trend could be observed for the travel time 
anomalies with epicentral distance. A unique model which would fit all the SH data for the 
region between south-eastern and north-western Europe could therefore not be obtained. 
Still, there are some features of the velocity structure valid for all paths. In the following I 
will show how these features of the upper mantle shear velocity structure are inferred from 
the data. 


SEISMOGRAMS WITH EPICENTRAL DISTANCES LESS THAN 20° 


The larger part of the data set consists of seismograms with an epicentral distance between 
15° and 20°. These data give the most detailed information of the upper mantle structure to 
the 400-km discontinuity. Interference of upper mantle phases can easily produce long and 
complicated wavetrains on which several arrivals may be recognized. For a structure without 
a low velocity zone (LVZ) the first arrival will be a refracted wave above the 400-km 
discontinuity. The second and third arrivals are then reflections from the 400- and 650-km 
discontinuity, respectively. This picture becomes more complicated when a LVZ is present. 
Then, the first arrival can be the Sn-phase, the phase travelling in the lid just above the LVZ. 
In the shadow of the LVZ, the Sn-phase is not present as a ray geometrical arrival, but may 
be observed as a diffracted phase. A later arrival is the phase which is refracted at the bottom 
of the LVZ, followed by the reflections of the two discontinuities. | mention only the two 
main discontinuities at this stage, because no evidence was found for other strong gradients 
or discontinuities. 

Some seismograms clearly present evidence for a high velocity lid above a LVZ. On the 
seismograms with an epicentral distance less than 17° (1900 km) a precursory arrival can be 
observed (see Fig. 5) which is identified as the Sn-phase. Its amplitude decays very rapidly 
beyond 17°, which would indicate the start of a shadow zone. However, the three seismo- 
grams with an epicentral distance less than 1900 km, showing the Sn-phase, are all recorded 
in the same station (NEO06). Unfortunately, there are no data from other stations in this 
distance range, but there is a seismogram from the same station at an epicentral distance of 
18.7° (2060 km) which also shows the Sn-phase. One can explain this observation in two 
ways. It is possible that the Sn-phase is not observed in the other stations because the 
shadow zone starts at a smaller epicentral distance along the ray-paths to these stations. This 
would imply a smaller lid thickness for the structure along these ray-paths than for the ray- 
paths to station NEO6. Another explanation is that the Sn-phase is not observed in the other 
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Figure 5. S-wave data on the transverse component of events in the eastern Mediterranean. The line is the 
predicted Sn arrival time (solid where observed, dashed where not). 


stations because of scattering. which diminishes its amplitude. Both interpretations, however, 
are indicative of heterogeneity to depths of at least 120 km. This is the depth of the start of 
the LVZ determined for the observations of station NEO6 for an average lid-velocity of 
4.4kms_ '. lt is not possible to determine the velocity structure in the lid accurately because 
of the sparseness of the data and the variations that these data already show, but, to be 
compatible with these and following observations, average velocities must lie within the 
range 4.44.7 kms !. The maximum thickness of the lid can vary between 70 and 110 km, 
depending on the lid-velocities assumed. 

Apart from the observations of the Sn-phase, the existence of a low velocity zone could 
also be established from arrivals from the underside of the low velocity zone, identified on 
seismograms with an epicentral distance larger than 21°. There are three events with such 
data. Fig. 6 shows the data for event 2-229, and the arrival from the gradient at the under- 
side of the low-velocity zone can clearly be identified. The other two events have source 

_depths of 25 and 61 km, respectively, and, because of the Jate arrival of the sS-phase, these 
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Figure 6. SH-wave data and travel-time curves of the S-phases of event 2-229. 


are less clear examples (Fig. 12a, b, 13a, b). These seismograms provide a very good 
constraint on velocity structure of the lower gradient of the LVZ, but the actual minimum 
velocity in the LVZ cannot be determined from body-wave data. 

The seismograms of event 3-186 show consistent waveshapes whose general features can 
be produced by a single model. The synthetics of model 1066B (Fig. 7a) for this event 
clearly show earlier arrivals from the 420-km discontinuity than observed in the data, while 
the reflection from the 670-km discontinuity is | or 2s earlier in the data. A model which 
includes a lid (max. depth 120km; v,= 4.5 kms 4), a LVZ (max. depth 300 km: min. 
v,=4.3kms’), about 4 per cent higher velocities than model 1066B just below the LVZ, 
and two discontinuities at 420 and 650 km depth explains the data reasonably well (Fig. 
7b). The velocity structure of the lid is not very well constrained. Its effect can be 
recognized by the first small negative onset, and it can be seen that this feature is not 
consistent between the data of NEO6 and the two other stations. The early arrival from the 
420-km discontinuity has been matched by an increase of the velocities below the LVZ by 
a constant value, but could also have been obtained by a decrease in depth of the 
discontinuity to about 360 km, or a combination of both effects. 

When the synthetics of event 4-042 (Fig. 8a) are calculated for the same model, a 
reasonable fit is obtained except for the data of NEO2. The amplitudes of most phases have 
not been modelled adequately, but the travel times generally agree within 2 s. A better fit to 
the data of NEO2 can be obtained by a small change in the velocity structure of the lid and 
upper part of the LVZ (Fig. 8b), but this model clearly does not fit the seismograms of the 
other stations so well, The first negative onset of the seismogram of NE06 is interpreted as 
the Sn-phase. It would require a decrease in the average velocity of the lid from 4.5 to 
4.4kms7' to model this phase better. There is some indication that the reflection from the 
650-km discontinuity is not matched adequately, but the phase is too weak or too 
contaminated by interference with other phases to justify any conclusive interpretation. The 
data at epicentral distances less than 20° do not constrain the structure below the 400-km 
discontinuity very well. 
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Figure 7. (a) SH-wave data of event 3-186 (solid line) and synthetics for model 1066B (dashed line): 
vertical lines represent the travel-time branches of the S-phases of the model. (b) Data and synthetics for 
model $3186. (c) Shear velocity models 1066B and $3186. 


BCO 


The phases of event 3-082, to the west of the two previous events, have their ray-paths 
for a large part beneath the Adriatic Sea. In contrast to the reasonable fit obtained for the 
seismograms of event 4-042, a bad fit is obtained for the data of event 3-082 using the same 
model that matches the data of event 3-186 (see Fig. 9a). Trial and error modelling leads to 
the conclusion that acceptable models require low velocities below the LVZ than the model 
for events 3-186 and 4-042. The subchannel velocities are then again in agreement with 
model 1066B in this depth range. The seismograms of NEO5 and NE04 (NE04 is not shown 
here because it is at the same epicentral distance as NEO5, and has a very similar waveshape) 
do not show an arrival from the bottom of the LVZ, or one coinciding with 400-km 
reflection (Fig. 9b). The large, broad, positive maximum on the seismogram of NEO6 can 
only be achieved when the travel-time branches of the 400-km reflection and the arrival 
from the LVZ are closely separated in time (Fig. 9b). Thus, the data of these three 
(including NE04) stations indicate that at an epicentral distance larger than or equal to 
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Figure 8. (a) SH-wave data of event 4-042 (solid line) and synthetics for model $3186 (dashed). (b) Data 
and synthetics for model $4042B. (c) Velocity models $3186 and S4042B. 


1950 km the phase from the 400-km discontinuity should arrive as a first arrival. However, 
the data of NEOQ3 can only be matched correctly when the phase from the LVZ arrives first 
(Fig. 9c), An earlier arrival of the phase from the LVZ can be achieved by an increase of the 
velocities in the lid and bottom side of the LVZ, while a late 400-km reflection, such as is 
observed, can either be explained by an increase in the depth of the discontinuity, or a 
decrease in the velocities below the LVZ. Unfortunately, the trade-off between these two 
possibilities cannot be resolved by the data of one seismogram only. 

The data of event 3-017, whose hypocentre lies very close to that of event 3-082, are very 
well modelled by the velocity structure obtained for event 3-082 (Fig. 10), although a 
slightly higher average lid velocity (4.45 kms! instead of 4.40 kms ') is adopted which 
yields a better travel-time fit of the total signal. This particular event has been modelled 
using the fault plane solution instead of the centroid moment tensor, as the former gives a 
better sS/S amplitude ratio for the data of station NEO6. 
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Figure 9. (a) SH-wave data of event 3-082 (solid line) and synthetics for model $3186 (dashed). (b) Data 
and synthetics for model §3082D. (c) Data and synthetics for model $3082C. (d) Velocity models $3186, 
$3082C, and $3082D. 


The data of event 3-218 very clearly point to heterogeneity to depths of approximately 
400 km. The main features of the data of stations NEO4 and NEOS5 are reasonably well 
matched by a model which has practically the same velocity structure as the model which 
fits the data of station NEQ3 for event 3-082 (the model of Fig. 9c; see Fig. 11a). However, 
the synthetics of this model do not match the data of station NEO? at all. This seismogram 
requires a velocity model which is more similar to that which fits the data of event 3-186 
(Fig. 11b). 
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Figure 10, SH-wave data of event 3-017 (solid line) and synthetics for model $3082D with an average lid- 
velocity of 4.45 kms”. 


SEISMOGRAMS WITH EPICENTRAL DISTANCES LARGER THAN 20° 


Seismograms with an epicentral distance larger than 20° provide more information about the 
structure of the upper mantle beneath the 400-km discontinuity. Depending on the velocity 
structure and the epicentral distance, several of the following phases may be observed: 
refractions from the bottom of the LVZ, transition zone, and lower mantle, and reflections 
from the two discontinuities. Most upper-mantle shear wave velocity models show that the 
phase from the lower mantle is the first arrival for epicentral distances greater than 25°. 
Phases retlected trom the 650-km discontinuity can be present to a distance of 30°, and 
from the 400-km discontinuity to 27°, which are the extreme values obtained for previous 
upper mantle models. 

When synthetic seismograms of event 4-173 are calculated for the model of event 3-082 
(for stations NEO4, NEO5, and NEO6), it can be seen that this model does not fit the data 
very well (see Fig. 12a). The ray-paths of the phases to station NEO4 for these two events 
differ little in azimuth, but the seismogram of event 4-173 gives a better constraint to the 
velocity structure within the transition zone because of its larger epicentral distance. The 
records of NEO4 and NE15 show a better fit with the synthetics for a model which has 
nearly the same velocity structure to the 400-km discontinuity as the model for event 
3-082 (stations NEO4, NEO5, NEOQ6), but has higher velocities below the 400-km 
discontinuity (see Fig. 12b). This model had earlier been obtained by matching the travel 
times of event 2-229 (Fig. 6). The data of NEO1 and NEO2 (Scandinavia) are not equally 
well matched by this model, which indicates that the velocity structure along the ray-paths 
to these two stations differs noticeably from that of the two stations in the Netherlands 
(NE04, NE15). 

A smaller gradient in the transition zone than given by model 1066B due to higher 
velocities below the 400-km discontinuity are also supported by the data of event 5-270. 
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Fig. 13(a) shows the data and synthetics for a model that matches the data of event 3-218 
for station NEQ2. A model with the same velocity structure in the transition region as for 
the previous event improves the fit of the first arrival, while a decrease of the velocities in 
the lid by 5 per cent improves the fit of the later arriving wavetrain. The synthetics for this 
model are shown in Fig. 13(b). 


INTERPRETATION OF THE S-WAVE MODELLING RESULTS 


Although the previous observations are indicative of distinct lateral heterogeneity, all upper 
mantle shear wave models show the same features: a high velocity lid (4.4—4.7 kms"), a 
low-velocity zone with minimum velocities of 4.20—4.38 km s”’ extending from 120 to 
300 km, and two discontinuities at about 420 and 650 km depth with velocity jumps from 
4.80 to 5.20 kms” and 5.48 to 6.1 kms ’, respectively. However, there are practically no 
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Figure 12. (a) SH-wave data of event 4-173 (solid line) and synthetics for model $3082D (dashed). (b) 
Data and synthetics for model $2229A. (c) Velocity models §2229A and §3082D. 


phases which give a constraint to the structure below the 650-km discontinuity, so this part 
of the model is not well resolved. This average structure is in many respects compatible with 
previous shear wave models obtained for Europe. Model M7 (Nolet 1977), which is inferred 
from higher mode Rayleigh wave data of Scandinavia and the western European Platform, 
shows a high-velocity lid and a low-velocity zone which are similar but less distinct than 
most of the models of this study. The 400-km discontinuity of this model is located at a 
shallower depth (380 km). Dost (1986), using data of the NARS array, showed that higher 
mode Rayleigh wave data for the western European Platform require a more pronounced 
low-velocity zone than that of the M7 model. Dost’s model (WEPL1) is remarkably similar 
to the average structure of the models obtained in this study to a depth of 300 km. The 
resolution of the surface wave data is not high enough at greater depths to allow a detailed 
comparison with the models from this study, Mayer-Rosa & Mueller (1973) were the first to 
present both P- and S-wave velocity models of the upper mantle beneath Europe. They 
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Figure 13, (a) SH-wave data of event 5-270 (solid line) and synthetics for model $3186 (dashed). (b) Data 
and synthetics for model §5270A. (c) Velocity models §3186A and S5270A. 


recognized changes in the velocity structure beneath south-east and south-west Europe, but 
explained the variations in the travel times for the two regions by models which include 
three low-velocity zones to a depth of 300km. Although many of these travel-time 
variations could also be interpreted by heterogeneity on a smaller scale, the most 
pronounced low-velocity zone in their models, at 200 km, seems to be a real feature. 

In contrast to previous studies, this analysis of broad band data presents a rough estimate 
of the degree of lateral heterogeneity in the shear wave structure between the eastern 
Mediterranean and north-west Europe to a depth of about 400 km. Unfortunately, there are 
not many seismograms which give us information about the structure in and below the 
transition zone, but these records do not indicate as much heterogeneity as is observed for 
shallower depths. One might have expected severe heterogeneity below the Mediterranean 
region, because of its tectonic complexity (McKenzie 1972). However, the modelling results 
show that heterogeneity at the receiver part of the ray-paths is also resolved: it is frequently 
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observed that one model cannot explain the features of all seismograms for one event. 
Although one should keep in mind that these models can only represent the average 
structure between source and receiver, many of these variations can be correlated with the 
tectonic setting along the ray-paths or with results of other studies. For example, 2-4 per cent 
higher velocities in lid and LVZ are needed to explain the S-wave data of station NEO2. 
when compared with the velocity structure determined for the seismograms of the other 
stations on the western European Platform (see the models for the events 3-218 and 4-042). 
This can be interpreted as a shield-like structure below northern Denmark, which is part of 
the Baltic Shield. Another remarkable observation is that, except for station NEO3, the 
seismograms of event 3-082 and 3-017 in the lonian Sea which have ray-paths beneath the 
Adriatic Sea require lower velocities in and below the LVZ, than the seismograms of 4-042 
in Greece. This observation is consistent with tomographic inversion results for the 
Mediterranean (Spakman 1986b) which show low P-wave velocities in the depth range from 
110 to 420 km beneath the Adriatic Sea when compared to the velocity structure east of the 
Adriatic and lonian Sea. Nolet, Panza & Worte!l (1978) determined a model tor the shear 
wave velocity in the lithosphere of the Adriatic subplate from surface waves. Their model 
shows an increase in velocity from 4.27 to 4.80 kms ' between depths of 35 and 80km. The 
gradient is larger, but the average velocity of the lid is equivalent with the mode! which is 
determined in this study. The depth to the LVZ is not well resolved by the study of Nolet 
et al. Another striking feature in the results of this study is the low value of the velocities in 
and below the LVZ obtained for the seismograms of event 3-218. There are presently no 
other studies which support this observation, but we will see that the same feature can be 
recognized in the P-wave velocity structure. 

It should be emphasized that the discontinuities have been modelled by a step in velocity 
only because the thickness of the discontinuities cannot be determined by this analysis. The 
reason for this is that amplitude and frequency information cannot be used due to the 
simplifications which are made with respect to attenuation. Only when attenuation can be 
modelled with a high accuracy, is it possible to determine the thickness of the discontinuities 
by waveform analysis of the reflected phases in studies of this kind. Furthermore, due to the 
trade-off between the velocity structure above the discontinuity and the depth of the 
discontinuity, it is not possible to determine the depth of the discontinuities with an 
accuracy greater than about 10 km. The depth of the discontinuities could be determined 
more precisely if the velocity structure above the discontinuities were known accurately. 
but, because of the heterogeneity to the 400-km discontinuity which is recognized in this 
study, a higher accuracy in a t-D model would be meaningless. 


6 Determination of the P-wave velocity structure 


Much of the information given in the previous section concerning the S-wave observations 
can also be related to the P-wave observations. Most of the seismograms which are used for 
S-wave modelling are also used for P-wave modelling. The data set is only difterent for those 
seismograms for which triggering occurred on the S-waves (no P-phases), or for which (one 
of the) horizontal seismographs were unstable (no S-wave modelling). Vertical-component 
seismograms were used to model the P-, pP-, and sP-phases. Similar to the result obtained for 
S-wave modelling, it appeared to be impossible to find one model that would fit all P-wave 


data. The P-wave structure of model 1066B could not be used as a starting model, since 
large travel-time differences between data and synthetics are observed (delays up to 8s), 


especially for seismograms with a small epicentral distance. The P-wave velocity structure 
of 1066B has a constant gradient from the crust to the 420-km discontinuity with an average 
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velocity which is obviously too high for central Europe, particularly at shallow depths. Nor 
do reference earth model PREM (Dziewonski & Anderson 1981) or other regional models fit 
the data well enough to be representative of the average velocity structure between south- 
east and north-west Europe. 


SEISMOGRAMS WITH EPICENTRAL DISTANCES LESS THAN 20° 


The P-wave structure to a depth of 400 km can be inferred from seismograms with an 
epicentral distance less than 20°. For models with a LVZ Pn-phases may be observed to an 
epicentral distance which depends on the velocity structure and depth extend of the lid. The 
first arrival for larger epicentral distances is ‘ne phase from the bottom of the LVZ. For 
velocity distributions without a LVZ phases refracted above the 400-km discontinuity arrive 
first. Phases reflected from the 400- and 650-km discontinuity may be recognized as Jater 
arrivals on seismograms with an epicentral distance less than 20°. 
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Figure 15. (a) P-wave data of event 3-218 (solid line) and synthetics for model P3218; arrows indicate the 
arrival times of the phases; TR represents the phase from the transition zone. (b) Data and synthetics for 
model P3218C. (c) Velocity models P3218 and P3218C. 


The P-wave data of this study do not show clear Pn-phases similar to the Sn-phases except 
for the data of one event. This could be an indication for the absence of a P-wave LVZ 
beneath the larger part of the region. However, there is evidence for an increase in the 
velocity gradient at approximately 200 km depth, which can be inferred from phases coming 
from this gradient zone. Moreover, it is observed that travel times of phases which bottom 
above the 400-km discontinuity are relatively slow. 

Fig. 14 shows the data of event 3-186 together with the synthetics for the tectonic model 
T9 (Burdick 1981), and for the model which is thought representative for the continental 
structure of north-west Eurasia, model K8 (Given & Helmberger 1980). The two models 
differ mainly in the velocity structure to a depth of 250 km. The tectonic model predicts a 
correct arrival time of the first phase, while the continental model gives a better fit to the 
670-km reflection (diffracted). It is observed that changes in the velocity structure below 
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Figure 16. (a) P-wave data of event 4-042 (solid line) and synthetics for model P3082F with a lid velocity 
of 8.1 km s‘!. (b) Data and synthetics for model P3082E with a lid velocity of 8.1 km s~’. (c) Velocity 
models P3082E and P3082F with lid velocities of 8.1 km s7?. 


250 km are needed to match the complete waveform. A model with a LVZ at the same 
depth interval as for the S-wave velocity model gives a better, although not optimal, fit to 
the data. 

Evidence for a P-wave LVZ comes from a precursory arrival on the seismograms of event 
3-218. This arrival has a high amplitude on the seismogram of NEO2 for this event, and is 
much smaller on the data of the other two stations. The later arrivals of the seismograms can 
very well be explained by phases which bottom at depths greater than 200 km, assuming a 
pronounced LVZ (see Fig. 15a, b). However, the first arrival can only be explained as a 
phase which has travelled through a high-velocity lid, or, equivalently, the Pn-phase. The 
WKBJ method cannot model this whispering gallery phase with the correct amplitude and 
frequency content without including a lot of internally reflected phases. The velocity 
structure that matches the data of NEOS does not fit the data of NEO4 equally well (only 
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Figure 17. (a) P-wave data of event 3-082 (solid Jine) and synthetics for model P3082E with a lid velocity 
of 8.1 km s~'. (b) Data and synthetics for model P3082F. (c) Velocity models P3082E with a lid velocity 
of 8.1! km s™' and P3082F. 


36 km difference in epicentral distance), and a bad fit is obtained for station NEO2 (see 
Fig. 15a). A velocity distribution with a slightly higher velocities above the 400-km 
discontinuity and a 20 km shallower discontinuity fits the data of NEO2 well. but does not 
match the data of station NEOS. Again it should be emphasized that in tnatching travel times 
there is a trade-off between the depth of a discontinuity and the velocity structure above the 
discontinuity. In this case it is clear that reflections from the 400-km discontinuity arrive 
earlier for a ray-path to a station in the north than for a ray-path to a more southerly 
station. Notice that the same result was obtained for the shear-wave velocity structure. The 
seismograms of this event require the most pronounced LVZ of the P-wave dataset, and this 
can be the reason why the Pn-phase is only observed for this event. 

P-wave data of event 4-042 are recorded in two stations in Scandinavia (NEO2 and NEO3), 
and one in France (NEO7). The P-wavetrain of the seismogram of NEO7 is matched 
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Figure 18, P-wave data of event 3-017 (solid line) and synthetics for model P3082F with a lid velocity of 
8.lkmst. 


adequately by a model that has a high-velocity lid (8.1 kms’), a LVZ (120-300km, 
min, ¥p = 7.85 kms~'), velocities below the LVZ which are 2.5 per cent lower than the 
model for event 3-186 (or models K8 and T9), and two discontinuities at 420 and 650 km 
(see Fig. 16a). This model does not give a good agreement between data and synthetics of 
the other two stations. A better fit to the seismogram of NEOQ3 is obtained by higher 
velocities below the lid (see Fig. 16b), but this velocity distribution still does not model the 
data of NEOQ2 adequately. The velocity structure below the LVZ of the model for NEO3 
corresponds to that of event 3-186. 

The P-wavetrain of event 3-082 recorded by station NEO3 is reasonably well matched by 
the model which was determined for event 4-042 for this station (Fig. 17a). The data of 
station NEO6, however, are better modelled by the velocity structure that was obtained for 
the seismogram of NEO7 for event 4-042 but with a slightly lower velocity lid (see Fig. 17b). 

The two seismograms of event 3-017, which has approximately the same hypocentre as 
event 3-082, are also well modelled by the velocity distribution of event 4-042 to station 
NEO7 (Fig. 18). As for the S-wave modelling, the fault plane solution has been used instead 
of the centroid moment tensor. The take-off angles of the P-phases are close to one of the 
nodal planes, which is not resolved accurately enough to model the first onset of station 
NEOS with the correct sign. The travel times, however, are satisfactorily matched, suggesting 
that the velocity model is adequate. 


SEISMOGRAMS WITH EPICENTRAE DISTANCES LARGER THAN 20° 


The seismograms with an epicentral distance greater than 20” have been used to investigate 
the structure below the 400-km discontinuity. The model which was determined for the data 
of station NEO3 for the events 3-082 and 4-042 gives a good fit to the seismograms of NE04 
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Figure 19. (a) P-wave data of event 4-173 (solid line) and synthetics for model P3082E. (b) Data and 
synthetics for model K8. (c) Velocity models P3082E and K8. 


and NE06 of events 4-173 (see Fig. 19a). The travel times of the seismograms of event 2-229 
of these two stations are also reasonably well matched by this model (see Fig. 20). The 
records of the stations NEO! and NEO2 of event 4-173, however, are clearly better modelled 
by model K8 (Fig. 19b). This could be an indication for higher average velocities in the 
upper 250 km of the mantle along the ray-paths to the stations NEO1 and NEOQ2. 

The model which gives a reasonable fit to most of the waveforms of the events 3-082, 
4-042, and 4-173 does not match the data of event 5-270 (see Fig. 21a). One way to improve 
the fit is to assume higher velocities in the lid and LVZ without changing the velocity 
distribution at greater depths (Fig. 21b). 

These results show that it is possible to match the P-wave arrival times by changes in the 
velocity structure above the 400-km discontinuity only. However, it should be emphasized 
that the resolution of the velocity structure below the 400-km discontinuity is not very high. 
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Figure 20, P-wave data of event 2-229 (solid line) and synthetics for model P3082E. 


INTERPRETATION OF THE P-WAVE MODELLING RESULTS 


The P-wave velocity structure for the region between the eastern Mediterranean and north- 
west Europe shows the same characteristics as the shear wave velocity distribution: a high- 
velocity lid (7.7 -8.3 kms7'), a low-velocity zone (120—300km) with minimum velocities 
varying between 7.55 and 7.90 kms~!. and two discontinuities at approximately 420 and 
650 km depth. The average velocity model below 300 km is similar to the models K8 (Given 
& Helmberger 1980) and T9 (Burdick 1981), which are thought representative for a shield- 
like and a tectonic province, respectively. The velocity structure of the lid and the minimum 
velocity in the low-velocity zone show some variations, but the average vaJues come close to 
those obtained by Mayer-Rosa & Mueller (1973): lid velocities ranging between 8.1 and 
8.2 kms™!, and a minimum velocity of 7.9 km s-' at 200 km depth. Model EWK for Eurasia, 
which is determined by England, Worthington & King (1977) trom travel-time and slowness 
measurements of earthquakes in the eastern Mediterranean and Adriatic regions, does not 
show a Jow-velocity zone. These authors note the low amplitudes of later arrivals between 
15° and 32° and observe much more scatter in the travel times of the phases when compared 
with data from Russian events, but these elements are not reflected in their model. 

This body wave analysis of broad band P-wave data presents evidence for heterogeneity in 
the P-wave velocity distribution of at least 2.5 per cent above the 400-km discontinuity. 
Several of the velocity variations which are inferred from the data are very well explained by 
the tectonic setting of the region which the phases cross. For example, the data of NEO and 
NEO2 are modelled by a velocity distribution with high velocities in lid and LVZ (event 
4-173). These high velocities are obviously related to the structure of the Baltic Shield (see 
e.g. Husebye & Hovland 1982). Another striking observation is that low velocities between 
the lid and 400-km discontinuity are required to model the phases of events 3-017, 3-082, 
and 4-042 which cross the Adriatic Sea (paths to NEO4, NEOS, NEQ6, and NEO7). The same 
observation was made for the shear wave structure and this feature is also recognized in a 
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Figure 21. (a) P-wave data of event 5-270 (solid line) and synthetics for model P3082E. (b) Data and 
synthetics for model P5270B. (c) Velocity models P3082E and P5270B. 


tomographic inversion by Spakman (1986b). Undoubtedly, the high velocities in the upper 
300 km which are obtained for the seismograms of event 5-270 are correlated with the cold 
subducted oceanic lithosphere beneath Crete. A profile of upper mantle P-wave velocity 
anomalies of this region shows this feature very clearly (Spakman 1986a, fig. 2). It is not 
possible to determine heterogeneity to depths greater than 400 km with this dataset due to 
its sparseness and the magnitude of the heterogeneity at shallower depths. 

As for the S-wave velocity models, the P-wave velocity models show first-order 
discontinuities. This is because amplitude and frequency information have not been used in 
the modelling to constrain the thickness of the discontinuity. Moreover, it is not possible to 
determine the depths of the discontinuities with high accuracy due to the trade-off between 
the velocity structure above the discontinuity and the depth of the discontinuity itself. 
However, depths of 420 and 650 km generally give the best fit to the data for representative 
values of the velocity distribution below the low-velocity zone. 
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8 Conclusion 


The upper mantle velocity structure beneath Europe has been investigated by modelling of 
body waves on broad-band seismograms. The data which have been used are P- and SH-waves 
of earthquakes in the eastern Mediterranean recorded by stations of the NARS array in 
north-west Europe. The seismograms have epicentral distances between 15° and 25°, and the 
phases sample the upper mantle to the 650-km discontinuity. [t is not possible to mode] the 
P- and S-wavetrains of the seismograms adequately by a single P- or S-wave velocity 
distribution. However, some features are consistent for the different P- and S-wave velocity 
models: a high-velocity lid, a LVZ and two discontinuities at depths of approximately 420 
and 650 km. The maximum depth of the lid is about 120 km and the velocities vary between 
4.4 and 4.7kms”’, and 7.7 and 8.3 kms! for P and S, respectively. Minimum values of the 
velocities within the low-velocity zone range from 4.20 to 4.38kms"', and 7.55 to 
7.90 km s~', respectively. The data show evidence for velocity variations of at least a few per 
cent in the depth range from 300 to 400 km, but do not require large velocity anomalies at 
greater depths. Although one might question the validity of each of the 1-D models, it is 
observed that many of the inferred velocity variations are related to the tectonics of the 
region between source and receiver. The depths of the discontinuities are not resolved very 
accurately (> 10 km), due to the trade-off between the depth of the discontinuity and the 
velocity structure above the discontinuity, which is clearly heterogeneous. Small scale 
undulations of the 650-km discontinuity, such as suggested by Paulssen (1985), can there- 
fore not be inferred by this kind of study. 

Broad-band seismograms are essential for a detailed study of the upper mantle by body 
wave analysis. Separate arrivals within the wavetrain of phases generally cannot be identified 
on long-period seismograms, which reduces the resolution which can be obtained for the 
upper mantle structure. Short-period seismograms lack the low-frequency information which 
is needed for S-wave modelling, and contain a high amount of high-frequency signals which 
are not related to the average upper mantle structure. The results of this study show that it is 
possible to resolve upper-mantle lateral heterogeneity on a regional scale by a detailed 
analysis of body wave phases. In principle, it would even be possible to apply tomographic 
techniques to the wavetrains of upper mantle phases. Unfortunately, the number of broad- 
band seismograms is presently too small to allow such a 3-D inversion. However, when 
combined with other information, the data can help to improve the resolution of the 3-D 
upper mantle structure. 
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Chapter 5 


Seismological evidence on the upper mantle 
discontinuities 


Direct information of the upper mantle discontinuities can be derived from 
seismological data. The seismic velocity and density contrasts across the transitions, as 
well as the depth and width of discontinuities are properties that are valuable to constrain 
models of the Earth’s upper mantle. 

In this chapter an overview is given of the seismological observations that have 
contributed to our present knowledge of the two upper mantle discontinuities (section 5.1). 
This chapter further presents the results of a study of upper mantle P-to-S converted waves 
that provides unambiguous evidence for a locally sharp 670-km discontinuity (section 5.2). 


5.1 AN OVERVIEW 


Seismological data that provide information about the upper mantle discontinuities can 
be classified in essentially two categories: seismograms with phases that sample the 
transitions by near-grazing incidence (regional data), and seismograms with waves that 
sample the transitions locally under near-vertical incidence. Section 5.1.1 presents a 
summary of the evidence inferred from regional P and S body wave studies, and section 
5.1.2 gives an outline of the information obtained from studies of upper mantle reflected 
and converted waves. 


Apart from a few minor additions sections 5.1.2 and 5.2 of this chapter are accepted for publication as: 


Paulssen, H., Evidence for a sharp 670-km discontinuity as inferred from P-to-S converted waves, J. Geophys. 
Res., in press, 1988. 
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5.1.1 Travel time and waveform analyses 


All regional velocity-depth models published in the two last decades, which are 
derived from travel time data or by waveform modelling, show two sharp velocity 
increases that are commonly referred to as the 400- and 670-km discontinuity. For 
simplicity, these transitions are generally displayed as first order discontinuities, but it 
should be emphasized that the data do not constrain these features with more resolution 
than 10 or 20 km (see e.g. Walck, 1984; Ingate et al., 1986). Variations in the depth of the 
discontinuities among these different (regional) models (see Nolet & Wortel, 1988 for a 
review) could reflect some genuine topography of these boundaries. On the other hand, 
there are substantial uncertainties in the depths of these transitions in the different velocity 
models. The velocity structure above the transition very must be known with a high 
accuracy in oreder to constrain the depth of the discontinuity with a high precision from 
travel times. As this is often not the case (it requires a large amount of data with phases 
that have their turning point just above the transition), it is likely that at least part of the 
depth variations among these models are artefacts of the interpretation. 

The inferred velocity (and density) contrasts across the transitions are of particular 
relevance for comparison with the values determined from high pressure and temperature 
mineral phases measured in the laboratory. The increase in compressional velocity across 
the 400-km and 670-km discontinuity is generally estimated to be in the range of 0.4 to 0.5 
km/s. There are fewer models of the upper mantle shear velocity, but these models indicate 
that the shear-velocity increase across the 400-km discontinuity is about 0.25 km/s, and that 
of the 670-km discontinuity is approximately 0.40 km/s. 


5.1.2 Reflected and converted phases 


Other evidence of the upper mantle discontinuities has been inferred from phases 
reflected or converted (from P to S, or S to P) at the transitions. 

Precursors of PKPPKP (P’P’) phases have been explained as underside reflections 
from the 400- or 670-km discontinuity, and they are generally denoted as P’400P’ or 
P’670P” (in this notation P’P’ would be P’OP’). A sharp 670-km transition has been 
suggested from observations of P’670P’ phases on short-period seismograms (Engdahl & 
Flinn, 1969; Whitcomb & Anderson, 1970; Adams, 1971; Bolt & Qamar, 1972; Teng & 
Tung, 1973; Fukao, 1977; Husebye et al., 1977, Nakanishi, 1986). Cleary (1981) 
questioned this interpretation and showed that such early precursors are also adequately 
explained by scattering of the PKKKP phase at the inside of the core-mantle boundary. 
More recently, Nakanishi (1986) argued that slowness estimates of the P’670P’ phases of 
several studies are not in agreement with those expected for PKKKP scattered phases. The 
small width of the transition has been inferred from the fact that short-period P’670P’ 
phases are indeed observed (assuming they are correctly interpreted). Richards (1972) 
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demonstrated that, to produce observable P’670P’ phases on short-period seismograms, the 
thickness of the transition must be less than 4 km (for a horizontally layered medium). This 
value has often been adopted in studies of this kind, and similar theoretical results were 
obtained by Teng & Tung (1973), and Fukao (1977). However, Richards (1972) also 
pointed to the importance of (even a slight) curvature of the reflector on the amplitude of 
the P’670P’ phase. In this regard, it is important to mention that P’670P’ phases are not 
consistently observed (Husebye et al., 1977; Nakanishi, 1986). Furthermore, P’670P’/P’P’ 
amplitude ratios vary among different studies (Teng & Tung, 1973). These observations 
can easily be explained in terms of focussing and defocussing effects of the discontinuity, 
although variations in the P’670P’/P’P’ amplitude ratio may also occur when one phase is 
hearer to a caustic than the other. In both cases it is likely that predominantly high- 
amplitude P’670P’ phases are observed, resulting in a overestimate of the reflectivity in the 
frequency band of short-period instruments, and, consequently, in an overestimate of the 
sharpness of the discontinuity (for acceptable values of the velocity contrast). 

Observations of P’400P’ phases on short-period seismograms have also been made, 
but they are relatively weak and appear less consistent than P’670P’ phases (e.g. Nakanishi, 
1988), which could indicate that the 400-km discontinuity is less sharp than the 670-km 
transition. 

First-order multiple reflections (’reverberations’) in the mantle of ScS phases (multiple 
ScS phases reflected once at one of the upper mantle discontinuities) have been used to 
determine the reflection coefficient (in the frequency band of long-period instruments) of 
the 670-km discontinuity for SH waves and the two-way (SH) travel time between the 400- 
and 670-km discontinuity beneath the western Pacific (Revenaugh & Jordan, 1987). 
However, these data do not constrain the fine structure of the discontinuities. 

Finally, studies of S-to-P and P-to-S converted phases have been carried out to 
investigate the nature of the discontinuities. Most of these studies have made use of long- 
period seismograms (P-to-S converted phases: Vinnik, 1977; Kosarev et al., 1984; Souriau, 
1986; Kind & Vinnik, 1988; S-to-P converted phases: Faber & Miiller, 1980, 1984; 
Baumgardt & Alexander, 1984) that not only demonstrated the presence of converted 
phases from the 400- and 670-km discontinuity, but also showed that large amplitude 
variations exist among these phases. These long-period data lack the high frequency 
information needed to obtain a high (spatial) resolution of the discontinuities, and only a 
few studies of short-period near source S-to-P converted phases have been made (Barley et 
al,, 1982; Bock & Ha, 1984; Richards & Wicks, 1987). These data suggest that the 670-km 
discontinuity is locally sharp beneath Tonga and Izu-Bonin, but also indicate that regional 
variations exist as systematic differences in the signal duration of the converted phases are 
recognized. However, these observations are not necessarily representative for an 
undisturbed upper mantle, since Tonga and Izu-Bonin are active regions of subduction, and 
the phases may even have converted at the subducted lithosphere itself. 

Summarizing the results of these studies, we must conclude that phases reflected or 
converted at the upper mantle discontinuities are often, but only intermittently observed. 
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The observations on short-period seismograms could be suggestive of a locally sharp 670- 
km discontinuity, but this inference has been disputed by Muirhead (1985) who argued that 
none of the interpretations is without ambiguity. The results of our study of P-to-S 
converted phases on short-period and broad-band data (section 5.2) gives to this date the 
most global and unambiguous evidence for a locally sharp 670-km discontinuity. From 
these data it furthermore appears that the 670-km discontinuity is more pronounced than 
the 400-km discontinuity. These results have recently been confirmed in a study of P-to-S 
converted phases beneath the Grafenberg array (Kind & Vinnik, 1988), a study of upper 
mantle P’dP’ phases by Nakanishi (1988), and observations of ScS mantle reverberations 
(Revenaugh & Jordan, 1987). 


5.2 EVIDENCE FOR A SHARP 670-KM DISCONTINUITY 


This section presents the results of a study of P-to-S converted waves from the two 
upper mantle discontinuities. Several aspects are of relevance to such a study and some of 
this background’ is presented in the following. 

It is well-known that (for a laterally homogeneous medium) a P wave incident on an 
interface not only generates transmitted and reflected P waves, but also transmitted and 
reflected S waves (see fig. 5.1). 


Figure 5.1. Ray geometry and notation of reflected, refracted, and converted waves generated by 
P waves incident on an interface at depth d. P waves are indicated by solid lines, S waves are 
dashed. 
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The P wave reflection coefficient and the conversion coefficients for the transmitted and 
reflected S waves are small for the velocity and density contrasts of the upper mantle 
discontinuities. Nevertheless, observations of upper mantle reflected and converted phases 
have been reported (see section 5.1.2), The transmitted phase which is converted at depth d 
and denoted Pads (see fig. 5.1), may be identified unambiguously by its SV polarization and 
travel time. Actually, any S wave arriving before the direct S phase must have been 
converted from P to S somewhere along its raypath. The time delay t between the direct P 
wave and the (slower) converted phase can be used to estimate the location of conversion 
(for a given velocity model). Furthermore, the P wave and the converted Pds wave arriving 
at the same location at the surface will have slightly different slownesses (see fig. 5.1). The 
values of the travel time difference t, the slowness difference Ap , and the ratio between the 
horizontal amplitude of the converted phase and the vertical amplitude of the P wave at 
72°, are given in table 5.1 for conversions at the 400- and 670-km discontinuity (for 
reference Earth models PREM, Dziewonski & Anderson 1981, and 1066B, Gilbert & 
Dziewonski 1975). 


Table 5.1: P-to-S converted phases (A=72°) 


400-km discontinuity 
model t(s) Ap (s/km) A‘,,/AP4s (%) 


rad! vert 
PREM 42 -0.0007 2.3 
1066B 44 -0,.0005 2.7 


670-km discontinuity 

model _(s)__Ap (skm) — AMg/Averr * (%) 
PREM 68 -0.0013 5.2 
1066B_ 69 -0.0013 8.2 


Table 5.1 shows that the amplitudes of the upper mantle converted phases are small. 
Identification of the conversions is further hampered by the signal related ’noise’ (the 
*coda’) of the P phase. It can be expected that stacking of the data for the appropriate 
slowness may enhance the signal-to-noise ratio of the Pds phase, whereas less coherent 
*coda noise’ is suppressed. Such a technique (more accurately described in section 5.2.2) 
has therefore been applied to the data selected for this study. 


5.2.1 Introduction 


The 400- and 670-km discontinuity mark the upper and lower boundary of the 
transition zone between upper and lower mantle. It should be kept in mind that the actual 
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depths of the discontinuities, generally not very well constrained, may differ from the 
values suggested by their names. In spite of the large number of studies on the upper 
mantle velocity structure, there is still no consensus about the nature of these two 
discontinuities. They can be related to phase transformations in olivine-rich mineral 
assemblages as well as to a change in chemical composition. Accurate seismological data 
are required in order to discriminate between these two possibilities, and to constrain 
petrological models of the upper mantle. So far, the (depth) resolution inferred from 
refracted seismic waves in the 15° to 30° epicentral distance range has proven insufficient 
to distinguish between a (presumably gradual) phase transformation and a (sharp) change 
in chemical composition. Detailed information on the sharpness of the discontinuities, 
which is the most important parameter to discriminate between the two, can be obtained 
from upper mantle reflected and converted waves, because they sample the discontinuities 
locally with relatively short wavelengths. A serious drawback of these phases is that they 
are generally of low amplitude. 

In this study the upper mantle discontinuities are investigated with near-station P-to-S 
converted waves from teleseismic events. These phases arrive in the coda of the direct P 
wave and have the polarization of teleseismic SV waves. Since the travel time of these 
waves is much less than that of the direct S wave, any SV-polarized wave in the coda must 
have been converted. Thus, contrary to other waves that have been studied, such as S-to-P 
conversions or underside reflections like P’670P’, their identification is not in doubt. The 
locality of conversion (the depth of the discontinuity) is constrained by the time delay with 
respect to P. In this study, the signal-to-noise ratio of the low amplitude converted phases 
is increased by stacking the seismograms of different events for a single station, which not 
only enhances the near-receiver phases ’generated’ near the station, but also allows an 
estimate of the coherence of the signal. 

It is shown that coherent energy, originating from the upper mantle beneath the 
Station, arrives in the time interval of 40 to 80 s after the P phase. Arrivals with an SV 
polarization, and a delay and slowness consistent with an interpretation of a P-to-S 
converted wave from the 670-km discontinuity (P670s phase), can be identified for several 
of the stations used in this study. These include five of the broad-band stations of the 
NARS array, the five short-period RSTN (Regional Seismic Test Network) stations in the 
U.S., and the short-period SRO station CHTO (Thailand). The results of this study not only 
confirm the results of previous investigations (Vinnik, 1977; Kosarev et al., 1984), but also 
complement these with high frequency information, presenting evidence for a locally sharp 
670-km discontinuity. The P670s arrivals observed for the stations of the NARS array are 
variable in signature, which may reflect a complex nature of the 670-km transition or a high 
degree of lateral heterogeneity beneath the stations. Upper mantle P-to-S converted waves 
from the 400-km discontinuity (P400s) are less consistently observed, and the data of this 
Study suggest that this transition is less sharp or of smaller velocity contrast than the 670- 
km discontinuity. 
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5.2.2 Data selection and method of analysis 


The data which have been used for this study are obtained from GDSN event tapes 
covering the period January 1985 - May 1986, and from NARS event tapes of the interval 
January 1983 - April 1987. The data have been selected according to the following criteria: 

1, Only short-period and broad-band data are considered. Long-period seismograms are 
discarded because of the inherent limitations in attainable resolution. 

2. No primary phases must arrive in the time interval of interest. This interval starts at 
about 40 s after the direct P wave when the amplitude of the crustal phases is negligible. 
To avoid interference with phases such as PP and PcP only teleseismic events are analysed 
(65° < A < 90°), and to prevent interference with pP and sP arrivals, only shallow (< 80 
km) and deep (> 400 km) events have been used. 

3. The available length of the seismograms must extend to a least 80 s after the P arrival to 
allow a possible identification of P-to-S converted phases from the 670-km discontinuity 
(P670s). 

4. The data must have a high signal-to-noise ratio and a clear P onset, which are checked 
by visual inspection. 

Furthermore, the number of seismograms per station must be sufficient to guarantee a 
reliable interpretation of the near-receiver phases of the P wave coda. This is the case for 
the short-period data of the 5 RSTN stations and SRO station CHTO (Thailand) with 15-47 
seismograms per station, and for the broad-band data of 5 NARS stations with 17-48 
seismograms per station (stations NE04, NEOS, NE06, NE15, and NE16). 

To allow an unambiguous interpretation of the near-receiver phases the data have been 
stacked. The stacking technique is essentially a t-Ap-stack, where tT in this case is the 
differential time with respect to the P phase and Ap the slowness difference with the direct 
P wave. The three components (vertical, radial, and transverse) of each seismogram are 
crosscorrelated with the waveform of the P phase on the vertical. The data are then 
normalized with respect to the autocorrelation of the P phase to give each seismogram 
equal weight. After that, the data are stacked using reference values of t and Ap 
appropriate for a P wave arriving at an epicentral distance of 72° for the data of the RSTN 
stations and station CHTO, and at 80° for the NARS data. These two values are 
approximately the mean values of the distribution of epicentral distances per station, and 
yield an optimum slowness resolution. The mean and standard deviation (o) of the stack 
(as a function of time) are calculated to obtain an estimate of the coherence in the stacked 
signal. Assuming a Gaussian distribution, the 95% confidence level of the signal (mean + 
20) can then be determined, allowing for an easy interpretation of the statistical 
significance of the stacked data. As an example, the t-Ap -diagram of the radial component 
of station RSCP with an indication of the 95% confidence interval (reflected in the width of 
the trace) is shown in fig. 5.2. Similar diagrams have been used to analyse the vertical, 
radial, and transverse component of the P wave coda of the stations. 

Phases for which the mean stacked amplitude is larger than twice the standard 
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Figure 5.2. t-Ap -stack of the radial component of the (correlated and normalized) seismograms 
of station RSCP (for an epicentral distance of 72°). Indicated by the width of the trace is the 
mean + 2 times the standard deviation (95% confidence interval) for each differential slowness 
relative to that of the direct P wave. The arrow points to a clear SV-arrival, and the star to the 
theoretical arrival of the P670s phase according to model 1066B. No primary phases arrive in the 
interval between 40 and 80 s the beginning of which is marked by the vertical line. 


deviation at that time are interpreted as ’significant’. In this respect, it is important to know 
that the standard deviation of the short-period and broad-band stacks generally is 1.5-3% of 
the amplitude of the P phase on the vertical in the 40-80 s interval. Arrivals with an 
amplitude smaller than about 3-6% of the direct P wave are therefore not expected to be 
resolved using the 95% confidence criterion. Phases have been interpreted as P-to-S 
conversions when ’arrivals’ on the radial component are not in phase with the signal on the 
vertical component, i.e. do not have a longitudinal polarization. 

The technique enhances all signals, on each of the three components, that are similar 
to the waveform of the P wave, and may therefore be used to investigate all phases that 
correlate strongly with the P-pulse. In most cases, the P-pulse has been defined by the first 
*sine-like’ waveform, usually about 2-3 s long. The method has been found to be quite 
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robust with respect to the selected signal provided that its duration is not too long (about 5- 
6s). This is probably due to the fact that crustal phases contaminate the P-waveform. The 
technique differs in this respect from the ’spectral ratio method’ often used to study the 
crustal structure beneath a station (e.g. Phinney, 1964; Langston, 1979; Owens et al. 1987). 
In the spectral ratio method, ’receiver functions’ for steeply incident P waves are obtained 
by the deconvolution of the radial component with the vertical component. It is then 
implicitly assumed that the vertical component contains the direct arrival (including source 
effects), but also most of the unwanted receiver crustal reverberations. The receiver 
function is therefore most sensitive to the crustal P-to-S converted phases, whereas 
reverberations are suppressed. This is an adequate technique for determining the crustal 
structure, but it is less satisfactory for the identification of (deep) upper mantle converted 
phases, because the incoming S (converted) phases are not affected in the same way by the 
crustal structure as the P wave. Hence, crustal multiples should preferably not be included 
in the ’source term’ as they degrade the deconvolution for the upper mantle converted 
phase. Using the technique described above, crustal multiples are enhanced instead of 
suppressed because of their similarity with the P-pulse. This explains the large amplitude 
of the ’ringing’ after the first P arrival as seen in fig. 5.2 for Ap = O km/s. (This 
information could be useful to further constrain the crustal P-wave models inferred from 
Station receiver functions that are most sensitive to the S-velocity structure.) 


5.2.3 Observations of P-to-S converted phases from the 670-km discontinuity 


Several ’arrivals’ of coherent energy can be recognized on the stack of fig. 5.2, but 
most interesting is the wavelet at T=68 s, Ap=-0.0025 s/km that is the only prominent 
arrival with the character of a near-vertical incident SV-wave in the time interval of interest 
(see fig. 5.3(a)). The phase at 63 s for a differential slowness of -0.035 s/km is clearly not 
an SV-arrival (see fig. 5.3(c)), whereas the phase at 71 s for the same differential slowness 
might be interpreted as a P-to-S conversion, although its polarization is more consistent 
with that of a P wave. The signal at t = 68 s and Ap = -.0025 s/km has a duration equal to 
that of the stacked P phase (at Ap=0 s/km) and its amplitude is 15% of that of the P wave 
on the vertical. To investigate the significance of the phase in more detail, the data set of 
this station is divided into two sets with different backazimuth (azimuth from the station). 
Although more noisy, both subsets show the same SV-arrival (see fig. 5.3(b)), indicating 
that there is no (or little) azimuthal dependence of the signal. 

Timing, amplitude, and polarization eliminate the possibility of a crustal origin of the 
phase. Furthermore, the arrival has a slowness smaller than that of the direct P wave, 
making an explanation as a crustal multiple more unlikely, as multiples have larger 
slownesses than the direct phase in a radially symmetric Earth. The polarization (near- 
vertical incident SV), and timing (coda of P wave) can only be explained by a phase 
converted from P to SV in the upper mantle beneath the station. Reference Earth models 
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Figure 5.3. (a) 95% confidence interval of the stacked (correlated and normalized) seismograms 
of the vertical (Z), radial (R), and transverse (T) component of station RSCP for a differential 
slowness of -0.0025 s/km. 

(b) 95% confidence levels of the radial component of the stacked seismograms of station RSCP 
with a backazimuth of 164-191° (upper trace, 8 seismograms), with a backazimuth of 319-326° 
(middle trace, 7 seismograms), and of total dataset (lower trace). 

(c) 95% confidence interval of the stacked (correlated and normalized) seismograms of the 
vertical (Z), radial (R), and transverse (T) component of station RSCP for a differential slowness 
of -0.0035 s/km. 


such as 1066B and PREM predict a delay of 68-69 s for a P-to-S conversion from the 670- 
km discontinuity (see table 5.1), suggesting that the phase is generated by the 670-km 
discontinuity. Its theoretical slowness is approximately 0.0013 s/km smaller than that of 
the P phase according to these models, a value slightly lower than estimated from the t- 
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Ap -stack, but the arrival clearly has a negative differential slowness as is predicted for a 
P-to-S converted phase. The reference models do not explain the high amplitude (15% of 
direct P wave on the vertical) of the P670s phase as observed on the stack of station RSCP. 
Calculated by the WKBJ method (for perfectly elastic, radially symmetric media), PREM 
predicts an amplitude of 5.2% of that of the P phase. Model 1066B, with a relatively high 
velocity (10%) and density (7%) contrast across the 670-km discontinuity (see Nolet & 
Wortel, 1988, for an overview of upper mantle velocity models), gives an amplitude of 
8.2% of that of the P phase. Obviously, the cause of the observed high amplitude must, at 
least partially, be sought in the structure beneath the station. 

Conjectures that local effects contribute to a strong P670s phase at RSCP 
(Cumberland Plateau, Tennessee) are confirmed by inspection of the individual 
seismograms. Fig. 5.4 shows examples of seismograms with arrivals identified as P670s. 
Enlargements of the relevant portions of the data of fig. 5.4(b-d) can be found in fig. 5.5(a- 
c). The amplitudes of the phases on the (short-period or broad-band) radial component 
reach values as high as 45-65% of that of the P phase on the vertical, but are, except for the 
example of fig. 5.4(a), still at about noise level which makes them difficult to identify. A 
notable feature of these seismograms is that the energy on the radial component is 
generally nearly as high as that on the vertical, even for the teleseismic P waves. The cause 
of this high amplitude signal on the radial component is not yet clear, but it might be 
related to location of the station on a sedimentary layer (see Owens et al., 1984, for the 
crustal structure beneath RSCP; this model predicts a radial/vertical amplitude ratio of 
about 25% for a teleseismic P wave with A = 72°, which is clearly not consistent with the 
observations). Bouchon & Aki (1977) modeled the response of a sedimentary basin to a 
local earthquake for wavelengths which are close to the dimension of the geological 
structure and a strong amplification of the horizontal motions is their main result. Although 
the cause of the high amplitudes on the radial component is not yet fully understood, it is 
evident that the amplitudes of the P-to-S converted phases at RSCP may not simply be 
interpreted by an exceptionally high P-to-S conversion coefficient. 

Another noteworthy feature of the high amplitude P670s phases is the similarity 
between the waveforms of the converted and the non-converted wave as shown in fig. 
5.5(a-c). The resemblance of the P and P670s waveforms is striking, given the noise level 
of the data (judged by the differences in the vertical and radial component for the P phase). 
The similarity between the P and P670s phase is also made clear by comparison of the 
amplitude spectra of the phases of fig. 5.4(c) (seen in fig. 5.5(d)). The wavelet of the 
converted phase resembles that of the direct P phase at least up to frequencies of 2 Hz 
(station noise dominates at higher frequencies). 

The short-period data of the other RSTN stations and station CHTO were analysed in 
detail to investigate whether the P670s phase could also be observed for these stations. 
However, the stacked seismograms of these stations did not yield equally conclusive 
results. There is some evidence for the presence of a P670s phase on the stacks of the 
stations RSON (Red Lake, Ontario), RSNT (Yellowknife, Northwest Territory), and RSSD 
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Figure 54. Vertical (Z), radial (R), and transverse (T) component of the short-period 

seismograms of stations RSCP of event 

(a) 1985-09-10: Depth: 33 km, Back azimuth: 329°, Delta: 66.8°, 

(b) 1985-11-12: Depth: 10 km, Back azimuth: 191°, Delta: 72.4°, 

(c) 1986-05-02: Depth: 33 km, Back azimuth: 325°, Delta: 89.9°, and 

(d) of the broad-band seismogram of event 1985-07-07: Depth: 32 km, Back azimuth: 167.8°, 

Delta: 69.3°. 


(Black Hills, South Dakota), where the signal is just statistically significant within the 95% 
confidence level (see fig. 5.6(a-c)). At t = 70, 72, and 69 s the amplitudes of the stacked 
signal are 4.3, 4.5, and 5.0% of that of the P phase for RSON, RSNT, and RSSD, 
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Figure 5.5. Waveforms of the direct P phase on the vertical and radial component (upper and 
middle trace), and of the P670s phase on the radial component (lower trace, 68 s shifted in time) 
of (a) event 1985-11-12 (cf. fig. 5.4b), 

(b) event 1986-05-02 (cf. fig. 5.4c), 

(c) event 1985-07-07 (cf. fig. 5.4d). 
(d) Frequency content of the P waveform on radial component (solid line), and of the P670s 


phase on the radial (dashed) of event 1985-11-12. (Vertical scale is linear). 


respectively, values that are in accordance with a PREM velocity contrast at 670 km. The 
interpretation of the data at RSSD is complicated by coherent signal on the vertical 
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Figure 5.6 (a-e). Stacked (correlated and nounagea) seismograms of the vertical (Z), radial (R), 
and transverse (T) component of the stations RSON, RSNT, RSSD, RSNY, CHTO, respectively. 
The number of seismograms, and differential slowness used for these stacks are indicated at the 
top of each diagram. 
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Figure 5.7. As figure 5.5 for the broad-band NARS stations NE04, NE0S, NE15, NE06, and 
NE16. Arrows point to ’arrivals’ on the radial component according to the 95% confidence 
criterion. 


component just before and after the possible P670s arrival. Note also that there are high 
amplitude signals on the transverse component of this station, which is probably due to the 
pronounced lateral variations in the crustal structure beneath this station (Owens et al., 
1987). There is no conclusive evidence for the presence of a P-to-S converted phase from 
the 670-km discontinuity on the stacked seismograms of station RSNY (Adirondacks, New 
York) and CHTO (Chiang Mai, Thailand), as can be seen from fig. 5.6(d-e). 

Similar to the short-period data, the broad-band data of the NARS array give varying 
results regarding possible P670s arrivals. Fig. 5.7 presents diagrams of the stacked 
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Figure 5.8. 95% confidence intervals of the stacked signal on radial component for the NARS 
stations of this study. Station name, number of seismograms used for the stack, and differential 
slowness are indicated at the top of each trace. Arrows point to positive ‘arrivals’. 


seismograms of these stations, where the arrows point to "95% confidence arrivals’ on the 
radial component. The data show a high degree of variability in the arrival times of 
possible SV phases, and in their ’signal-to-noise ratio’ (the radial components have been 
drawn to the same scale in fig. 5.8). The most consistent arrival is observed for t = 66 s, 
while another less coherent SV arrival might be recognized at a delay of approximately 72 
s. In the following, we will therefore concentrate on these wavelets. 

The stacked seismograms of station NE15 (Valkenburg, The Netherlands) show a 
clear SV arrival at t = 66 s, with an estimated slowness of Ap = -0.0010 s/km, and an 
amplitude which is 7.3% of the amplitude of the direct P phase on the vertical with a 
standard deviation of 2.4% (7.3 + 2.4 %, this notation will be used in the following). The 
phase has a duration equal to that of the stacked P phase and is observed on different 
substacks of the data. It therefore seems to be a coherent phase, but it is difficult to identify 
on the individual seismograms. Apparently, the phase is generally not of high amplitude, 
except perhaps for the one shown in fig. 5.9(a) with a delay of 65 s. The other two stations 
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Figure 5.9 (a-c). Examples of P670s phases for different NARS stations in northwest Europe. 
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in the Netherlands (NEO04 Witteveen, NEOS Utrecht) and station NE06 (Dourbes, Belgium) 


show at best poorly resolved SV waves at 66 s with amplitudes of about 5.0 + 2.3 


%. 


Examination of subsets of the data reveals that the arrival at 66 s is not consistently 
observed at NEO4 and NEOS (see also Paulssen, 1985), but, although with smaller 
amplitude, appears more coherent at NEO6. Inspection of the individual seismograms 
shows that clear arrivals with a consistent delay of approximately 66 s can be identified on 


some of the records. Fig. 5.9 shows examples of phases that have been interpreted as P- 


to- 


S conversions from the 670-km discontinuity for different NARS stations in northwest 
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Figure 5.10. (a) Example of a P670s phase for station NE04. 


(b) Waveform of the P phase on the vertical component (solid line) and the P670s phase on the 
radial component shifted by 66 s (dashed). 


(c) Amplitude spectra of the P and P670s waveforms. 
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Figure 5.11. Stacked seismograms of NE04 and NEOS for Ap = 0 s/km showing an SV arrival at 
T=80s. 


Europe. Fig. 5.10(a) has been taken from Paulssen (1985) to show the similarity (fig. 
5.10(b)) for a long and high-frequency source time function. The converted and non- 
converted phase resemble each other up to a frequency of about 1 Hz, the cutoff-frequency 
of the NARS instrument response (see fig. 5.10(c)). 

The most coherent signal with an SV character observed for station NE16 (Clermont 
Ferrand, France) is at the arrival at 72 s (Ap = -0.0010 s/km, ampl.= 8.7 + 2.7 %). 
However, the phase is not recognized on any of the individual seismograms due to the 
generally high amplitude (’noisy’) coda of the data. A similar remark could be made for 
the arrival at 73 s for station NEO4, although the ’coda noise’ is of this station is smaller in 
amplitude. The wavelet at 70-76 s observed for station NE06 is a feature of seismograms 
with a northeast backazimuth only. Although some of the seismograms of these stations 
show a high amplitude signal on the radial component at about 72-74 s after the P arrival, 
these cannot unambiguously be interpreted as upper mantle converted phases. 

Another interesting feature, shown in fig. 5.11, is a very clear SV arrival at 80 s in the 
coda of NE04 for data with a backazimuth of about 345° to 41° (Ap = 0 s/km, ampl.= 11.4 
+ 3.1 %), and a similar but less clear phase for NEOS (Ap = 0 s/km, ampl.= 10.1 + 4.6 %) 
for which the available seismograms have a predominantly northeast backazimuth. The 
phase is not observed for any of the other stations, and is recognized on the seismograms of 
NE04 and NEOS by an increase in amplitude. Waveform comparisons with the direct P are 
not as convincing as for phases arriving at approximately 66 s, the phase generally being 
recognized by an increase of energy on the radial component. This, together with the 
estimated differential slowness of 0.0 s/km, and the limited number of stations and 
azimuthal range for which these observations are made, hampers the interpretation as 
simple upper mantle P-to-S converted waves. 


84 


RSNT, N=37, dp=-.0020 s/km 


47 48 49 50 51 
Time (s) 


Chapter 5 


RSNY, N=37, dp=-.0005 s/km 


46 47 48 49 50 
Time (s) 


NE04, N=48, dp=-.0040 s/km NE06, N=37, dp=-.0025 s/km NE16, N=38, dp=-.0010 s/km 


40 44 48 40 44 48 52 40 44 48 52 
Time (s) Time (s) Time (s) 


Figure 5.12. Stacked seismograms with possible P-to-S conversions from the 400-km 
discontinuity. 


5.2.4 Observations of P-to-S converted phases from the 400-km discontinuity 


The observations of the P670s phase stimulated a detailed investigation into the 
presence of other converted phases from the upper mantle. Although there is some 
evidence for the existence of other SV arrivals in the time window of 40 to 80 s after the P 
phase, none of these phases is very systematically observed. Most coherent and clear are 
the phases with a delay of 43-49 s. Two RSTN stations show an SV arrival at 48-49 s after 
the P phase (RSNY: t = 48 s, Ap = -0.0005 s/km, ampl.= 5.0 + 1.3 %; RSNT: tT = 49 s, Ap 
= -0,0020 s/km, ampl.= 4.5 + 1.2 %), and three NARS stations show an SV phase with a 
delay of 44-46 s (NE04: t = 44 s, Ap = -0.0040 s/km, ampl.= 6.3 + 1.9%; NE06: t = 46 s, 
Ap = -0.0025 s/km, ampl.= 2.8 + 1.4%; NE16: t = 44 s, Ap = -0.0005 s/km, ampl.= 9.0 + 
2.7%). Fig. 5.12 shows the 95% confidence intervals of these data. If these phases are 
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interpreted as P-to-S converted phases from the upper mantle they must be related to the 
400-km discontinuity. Model PREM and 1066B predict delays of 42 and 44 s for 
discontinuities at depths of 400 and 420 km, respectively, and delays of 43 and 46 s if a 
more acceptable continental crust is incorporated. The differential slowness is -0.0007 
s/km according to PREM and -0.0005 s/km for 1066B, and theoretical amplitudes are 
approximately 2.5% of the P wave. The amplitudes of the P400s phases are smaller than 
those of P670s phases, mainly because of a smaller velocity contrast in the two reference 
models at the 400-km discontinuity. Observations of P400s phases are indeed smaller in 
number than those of P670s phases, but the stacked amplitudes range from 2.8 to 9.0 %, 
values that are comparable to those of P670s phases. These variations in estimated 
amplitude, travel time, and slowness, point to lateral heterogeneity at or above the 400-km 
discontinuity. No clear P400s phases have been identified on individual seismograms, 
which could imply that there are no high-amplitude P400s phases, or that there is more 
waveform distortion than for the P670s phases. More data are needed to resolve this 
question, and to allow an interpretation in terms of the width and velocity contrast of the 
discontinuity. 


5.2.5 Interpretation 


In the previous section I showed that the stacked seismograms allow an identification 
of the upper mantle P-to-S converted phases beneath the seismic station if the signal-to- 
noise ratio is sufficiently increased by stacking. The amplitude estimates as obtained from 
the stacks should be interpreted as lower bounds of the ’true’ mean amplitudes, since small 
travel time differences due to crustal and upper mantle heterogeneity beneath the station 
will reduce the amplitude of the stacked signal, especially for the short-period data. 
Furthermore, it should be kept in mind that the inferred differential slownesses are only 
estimates of which the resolution depends on the distribution of epicentral distances. The 
visually determined ’slowness resolution’ of the t-Ap-stacks is for all stations 
approximately 0.0005 s/km, except for RSON (+ -0.0010 s/km). 

The results of this study show clear P-to-S converted phases from the 670-km 
discontinuity for the short-period data of station RSCP (t = 68 s, ampl.= 13.7 + 3.6%, Ap = 
-0.0025 s/km), and for the broad-band data of station NE15 (t = 66 s, amp1.= 7.3 + 2.4%, 
Ap = -0.0010 s/km) and NE16 (t = 72 s, ampl.= 8.7 + 2.7%, Ap = -0.0010 s/km). 

Less clear are the P670s phases on the short-period stacks of RSON (1 = 70 s, ampl.= 
4.3 + 1.7%, Ap = -0.0020 s/km), RSNT (t = 72 s, ampl.= 4.5 + 2.1%, Ap = -0.0010 s/km), 
RSSD (t = 69 s, ampl.= 5.14 2.5%, Ap = -0.0005 s/km), and possibly the phase at RSCP (t 
= 71 s, ampl.= 11.0 + 4.9%, Ap = -0.0035 s/km). The other broad-band NARS stations 
present a picture of signals with an SV wave character between 66 and 76 s. These signals 
are just statistically significant within the 95% confidence limit, but it is difficult to trace 
their origin. The arrivals can be interpreted as conversions from interfaces at different 
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depths, but they can also be explained by lateral heterogeneity beneath the station, resulting 
in variations of the arrival times for phases generated by the 670-km discontinuity. The 
data of station NE06 seem to point in this direction as the wavelet between 70 and 76s is a 
feature of data with a northeast backazimuth only. However, most clear P670s arrivals 
have been identified on seismograms of different NARS stations in northwest Europe with 
a delay of 66 + 2 s, and this arrival is more or less consistently observed on the stacks of 
these stations (NE04, NEOS5, NEI5, and NEO6) for different azimuths. The question of 
whether the later arrivals are due to different interfaces or to lateral heterogeneity cannot be 
resolved with this data set. The intermittent nature of observations of P-to-S converted 
phases from the 670-km discontinuity is further inferred from the short-period stacks of 
station RSNY and CHTO which do not show any evidence for such a phase exceeding the 
standard deviation of 1.3% and 1.9% respectively. 

The identified P670s phases on the short-period and broad-band seismograms of 
station RSCP and on the broad-band seismograms of NARS stations show extremely high 
amplitudes that cannot be explained by (acceptable) radially symmetric upper mantle 
models. It is believed that the high amplitudes are at least partly caused by amplification 
effects in the local structure beneath the station. Supporting this view is the observation 
that clear P670s phases on individual seismograms are predominantly identified for stations 
located on sediments, which might imply that focussing effects of the sediment-basement 
interface play an important role. However, high-amplitude P670s phases can also result 
from focussing effects of other structures in the crust or upper mantle. In particular, 
focussing by warping of the 670-km discontinuity is very effective due to its large 
impedance contrast for P-to-S converted waves. The travel time and slowness variations as 
observed for the P670s phases are further suggestive for topography of the 670-km 
discontinuity, although effects of upper mantle heterogeneity may not be neglected. Since 
a large degree of upper mantle lateral heterogeneity is likely to be present beneath different 
regions (particularly for western Europe: Spakman et al., 1988; see also chapter 4), it is not 
appropriate to interpret the observed travel times, slownesses, and amplitudes in a 
quantitative way. Nevertheless, to get an idea of the magnitudes of the parameters 
involved, it is simply calculated that for a lens-shaped topography of the 670-km 
discontinuity a radius of curvature of about 1300 km is needed to double the amplitude of 
the converted phase (in an otherwise homogeneous upper mantle). This value corresponds 
to a change in depth of the transition of 15 km over a horizontal distance of 200 km, or, 
equivalently, to a travel time variation of about 1.5 s for the P670s phase. 

More important for an interpretation of the nature of the 670-km discontinuity is the 
waveform similarity of the converted phases to the direct P wave. The waveforms of the 
converted phases resemble those of the P wave signal up to 1 Hz for the broad-band data 
(see figures 5.4(d), 5.9, and 5.10), and for the short-period data in one particular case even 
up to 2 Hz (fig. 5.5(d)). This implies that the 670-km discontinuity is, at least locally, 
extremely sharp. The transfer functions shown in fig. 5.13 display the frequency-dependent 
amplitude behaviour of converted phases for different widths of the discontinuity. They 
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Figure 5.13. Transfer functions, of the frequency dependence of P670s phases for different 
widths of the discontinuity. The curves represent the (frequency dependent) amplitude ratio of 
phases converted at a transition of finite width compared to those converted at a first order 
discontinuity. Indicated for each curve is the width of the transition in km. 


are calculated by the propagator matrix method, assuming a ‘linear gradient’ within the 
transition which is approximated by a sequence of layers with a thickness of 1 km each. 
Assuming that the high-frequency content has decayed by 40 to 50% relative to the low 
frequencies, it can be inferred that the width must be less than 5 to 6 km for waveform 
*similarity’ up to 1 Hz (beneath station NE04, fig. 5.10(c)), or less than 2.5 to 3 km for 
waveform ’similarity’ up to 2 Hz (beneath station RSCP, fig. 5.5(d)). This observation has 
important implications for petrological and geodynamical models of the mantle (see 
discussion and chapter 6). 

Observations of P-to-S converted waves from the 400-km discontinuity on the stacked 
seismograms are smaller in number, and no clear P400s phases have been identified on 
individual seismograms. This may indicate that there are no high amplitude P400s phases, 
or that there is more waveform distortion for the P400s phases than for the P670s phases. 
This last interpretation could imply that the 400-km discontinuity is less sharp than the 
670-km discontinuity. However, to produce an amplitude of 4-5 % of the P wave ona 
stack of short-period data, the 400-km discontinuity must be both sufficiently sharp (< 10 
km) and define a reasonable velocity contrast (> 5% if no focussing effects are included). 
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The most important result of this study is the evidence for a locally sharp 670-km 
discontinuity from identifications of P670s phases on seismograms. A sharp 670-km 
discontinuity has previously been suggested from observations of precursors to P’P’ phases 
and near source S-to-P converted phases, but these results were not without ambiguity (see 
section 5.1.2). In this respect, the waveform comparisons of P and P-to-S converted phases 
(Paulssen, 1985; this study) leave no doubt. The number of observations of upper mantle 
reflected and converted phases, and the estimated locations of reflection or conversion 
indicate that the discontinuity is a world wide feature present under continents and oceans, 
subduction zones and ridges. There is growing evidence that the transition is locally sharp 

These seismological observations have important geodynamical implications. Lees at 
al. (1983) have shown that, for plausible mantle mineralogies, models that invoke phase 
transitions alone produce a reflectivity in the short-period frequency band that is an order of 
magnitude smaller than inferred from P’670P’. The reflectivity for short periods critically 
depends on the width of the phase transformation. Using the high pressure and temperature 
data compiled by Jeanloz & Thompson (1983), Lees et al. calculate that phase transitions 
occur over a 10 to 40 km depth interval, depending on the mineral assemblage. The values 
in this width range are too high to explain the seismic observations. Recent experimental 
work by Ito & Takahashi (1987) shows that the dissociation of y-spinel into perovskite plus 
magnesiowustite, probably the most important phase transition at a depth of 650 km, is 
completed within a pressure interval of 1 GPa (= 25 km). Although this is sharp in a 
petrological sense (for a divariant phase transformation), it should be resolved with an 
accuracy of less than 0.2 GPa (5 km) to be able to compare the petrological results with the 
seismic observations. The dissociation of majorite, a garnet solid solution, evidently cannot 
explain the seismic 670-km discontinuity, because it occurs over a relatively wide pressure 
interval of about 4 GPa (Ito & Takahashi, 1987). A chemical transition, or a combination 
of a chemical and phase transition, can easily explain the sharpness of the 670-km 
discontinuity (see also Lees et al., 1983). A compositional change across the discontinuity 
would most naturally imply a chemically distinct upper and lower mantle, possibly 
involving 2-layered convection (Richter & Johnson, 1974; Christensen & Yuen, 1984). An 
alternative model is proposed by Ringwood & Irifune (1988). In this model, the upper and 
lower mantle are assumed to be of similar bulk composition separated by a chemical 
boundary layer at a depth of 600 to 700 km, which is supposed to have formed by effects of 
phase changes and differentiation in subducted slabs. The phase transformation of y-spinel 
to perovskite plus magnesiowiustite together with the chemical change of (lenses of) former 
basaltic crust on top of ancient oceanic lithosphere could account for the sharp 670-km 
discontinuity and the variability of the (seismic) observations. 

At present it is not clear whether the upper and lower mantle are chemically distinct or 
not. More seismological studies are required to improve our knowledge of seismological 
aspects of the 670-km discontinuity to be able to answer some basic questions: Is the 670- 
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km discontinuity only locally sharp, or is its sharpness a worldwide feature? Is topography 
of the 670-km discontinuity, as suggested by Richards & Wicks (1987), only locally 
significant or also on a global scale? What is the origin of the multiple SV-arrivals in the 
66-72 s interval (this study)? What is the exact’ velocity and density contrast across the 
670-km discontinuity? 

Short-period and broad-band body wave studies of converted and reflected phases now 
provide the most detailed information about the sharpness and depth of the 670-km 
discontinuity, but they could be complemented with those of long-period (or low-passed 
broad-band) data to investigate the frequency dependence of the variations among the 
observations. On the other hand, upper mantle refracted phases are more adequate to 
determine the velocity contrast across the transition. Although we are still far from ’ global 
mapping’ of the upper mantle discontinuities, an increased station density combined with a 
sufficient quantity of broad-band digital data would contribute to solve many of the 
ambiguities that still exist in detailed upper mantle studies. 

Apart from seismological information, major advances in our understanding of the 
670-km discontinuity are expected from other fields, especially from experimental 
petrological studies under high pressure and temperature conditions (see chapter 6). 
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Chapter 6 


The upper mantle 


Not only seismology has contributed to present insights into the constitution of the 
Earth’s upper mantle. Other geophysical data, as well as geochemical and geological 
observations, have provided valuable information about the structure of the upper mantle. 
The observations are widely varying in character, yielding information about different 
aspects of the upper mantle. In this chapter, a summary is given of some of this evidence 
and several of the models that have been proposed for the upper mantle. It is not meant to 
present a ’complete’ picture of our knowledge of the upper mantle, but it intends to give 
some insight into phenomena that are well established and controversies that to this date 
still exist. 


6.1 PETROLOGY OF THE UPPER MANTLE 


Petrological models can be seen as essential concepts of the mantle because they have 
implicit bearings on thermal and geodynamical aspects of the mantle and on aspects of the 
Earth’s early evolution. Thus, a uniform mantle composition would suggest no chemical 
differentiation after, or concurrent with, the Earth’s accretion, and would permit a regime 
of mass transport that encompasses the entire mantle. Petrological data from laboratory 
experiments have been used to investigate the viability of compositional models of the 
mantle and several aspects of this discussion will be presented in section 6.1.3. But first, in 
sections 6.1.1 and 6.1.2, a brief overview of the presently available petrological 
information of the upper mantle is given. 


6.1.1 Important phase transformations 


It is generally accepted that the most important upper mantle constituents are olivine 
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((Mg,Fe)2SiO,), pyroxene ((Mg,Fe)SiO3), and garnet ((Mg,Fe,Ca)3Al,Si,0,). Major 
phase transformations in these components occur in the 10-30 GPa pressure interval, and 
this interval corresponds well with the 300-800 km depth range over which the seismic 
discontinuities and gradients are observed. Since phase changes in these systems could 
possibly account for all density and velocity transitions in the mantle, they have extensively 
been investigated, particularly in the last 10 to 15 years. Thus, it is known that 
(Mo 9F€p,;)25iO,-olivine converts to the B-spinel structure at a pressure of about 13 GPa 
(400 km), which subsequently transforms to the y-spinel phase (at about 15 GPa, 440 km), 
and to an assemblage of (Mg,Fe)SiO3-perovskite and (Mg,Fe)O-magnesiowistite (at about 
23 GPa, or 650 km). (Mgo9Fep)SiO3-enstatite transforms to an assemblage of B- 
(Mg,Fe)2SiO, plus stishovite, to y-(Mg,Fe)2SiO, plus stishovite, to (Mg,Fe)SiO3-ilmenite, 
and to (Mg,Fe)SiO,-perovskite with increasing pressure. Most of the phase 
transformations in the pyroxene system occur at approximately the same pressures as in the 
olivine system, except that little is known of the stability field of ilmenite (for a review see 
Jeanloz & Thompson, 1983). In the presence of garnet, pyroxene will gradually dissolve in 
the garnet structure with increasing pressure, giving rise to a garnet-majorite solid solution. 
Irifune (1987) recently reexamined this pyroxene-garnet transition, because there was no 
general agreement on the details of this transformation. He showed that pyroxene 
gradually dissolves in the garnet structure up to pressures of 13-14 GPa, forming a single 
phase garnetite at 16 GPa (S70 km). Ito & Takahashi (1987) showed that majorite-garnet 
further dissociates into Ca- and Mg-rich) perovskites (and possibly garnet) over a broad 
(2-4 GPa) pressure interval starting at temperature and pressure conditions corresponding 
to a depth of about 600 km. Such a broad phase loop is in contrast to the transformation of 
y-(Mg,Fe).SiO, to perovskite which is completed over a pressure interval of less than 1 
GPa (Ito & Takahashi, 1987). 


6.1.2 Are the upper and lower mantle chemically distinct? 


Although compositional changes across the 670-km discontinuity had been proposed 
earlier, Liu (1979a,b) was the first to suggest this on the basis of results from high pressure 
and temperature laboratory experiments. He argued that the change in density and bulk 
sound velocity (vp-4/3v2)* across this transition must be the consequence of a change in 
composition from an olivine-rich upper mantle to a predominantly perovskitic lower 
mantle. However, uncertainties in the thermal expansion coefficient (a) and bulk modulus 
(K) of the perovskite phase affect the reliability of his interpretation (see e.g. Jackson, 
1983). Knittle et al. (1986) and Knittle & Jeanloz (1987) recognized the importance of 
well constrained physical properties of the perovskite phase for the petrological 
interpretation of the lower mantle, and carried out diamond-anvil laboratory experiments to 
determine a, K, and K’ (the pressure derivative of K) of (Mgo9Fep,;)SiO3-perovskite. They 
showed that the zero-pressure density of the lower mantle is adequately explained by a pure 
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perovskitic composition, whereas a pyrolitic material in a lower mantle mineralogy yields a 
(zero-pressure) density that is 2% smaller than that of the lower mantle. These results 
would be in support of chemical stratification of the mantle, for instance by iron or silica 
enrichment in the lower mantle. Although these data provide at present the most 
compelling petrological arguments in support of a chemically distinct upper and lower 
mantle, there may be uncertainties, especially in the extrapolation of the data to zero 
pressure. 

Lees et al. (1983) have addressed the subject from a different point of view. These 
authors compared the reflection properties obtained from the P’670P’/P’P’ amplitude ratio 
with those calculated for phase transition and compositional change models. Phase 
transitions in the olivine, pyroxene, and garnet systems are predicted to occur over a depth 
interval of 10 to 40 km (on the basis of the available data at that time). A transition over 
such an interval is too wide to explain the high amplitudes of the the observed P’670P” 
phases, or the observations of P-to-S converted phases presented in chapter 5. Lees et al. 
further showed that acceptable quantities of iron or silica enrichment across the 670-km 
phase transformation (or in the high-pressure phase regime) would produce a 
*discontinuous’ transition of which the calculated reflection coefficient is consistent with 
the seismic data. However, it should be noted that the amplitude of a seismic phase may be 
affected significantly by focussing effects along an individual raypath, and this may 
invalidate a simple petrological interpretation of the seismological amplitude information. 

In this regard, it is important to note that unambiguous evidence for a locally sharp 
670-km discontinuity is presented in chapter 5. It is clear that the observations of P-to-S 
converted phases are readily reconciled with a change in chemical composition. The 
important question still to be answered is whether a phase transition can occur over a depth 
interval of a few kilometers only. Recent results of Ito & Takahashi (1987) showed that the 
y-spinel to perovskite plus magnesiowistite phase transformation at high temperatures 
(1600° C) occurs over a small pressure interval, and this could alter the conclusions of Lees 
et al. (1983). They constrained the width of this phase loop to a pressure interval of less 
than 1 GPa. This is still 5 times larger than the thickness implied by the seismic data. At 
present, it seems unlikely that the transformation occurs over an interval of less than 0,2 
GPa (~ 5 km) (see also Jeanloz & Thompson, 1983). The majorite-perovskite phase 
transition must be excluded as a possible explanation of the sharp seismic 670-km 
discontinuity because it occurs over a wide pressure interval of about 2 to 4 GPa (Ito & 
Takahashi, 1987). 

The only way to reconcile a sharp 670-km discontinuity with a broader divariant phase 
transformation may be to invoke nonequilibrium effects caused by large vertical motion, as 
suggested by Loper (1985). Material moving fast through a phase boundary may then not 
immediately form new crystals at equilibrium conditions. ’Overstepping’ will occur when 
the pressure and temperature conditions are such that the threshold energy for the phase 
reaction in the dynamic situation is reached (Schuiling, pers comm., 1988). However, the 
importance of this effect for the y-spinel to perovskite plus magnesiowiistite phase 
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transition at mantle conditions is not yet known. We conclude that, the most 
straightforward interpretation of the seismic 670-km discontinuity presently is a transition 
of a (at least partially) compositional nature. 


6.1.3 Compositional models of the upper mantle 
PYROLITE MODEL 


Ringwood (1975) presented a review of the experimental petrological data at that time 
available, and related these to upper mantle petrological, mineralogical, geochemical data, 
and cosmochemical and seismic evidence. Ringwood suggested that the main constituent 
of both upper and lower mantle is a theoretical mineral assemblage called ’pyrolite’. It has 
the composition of what is thought to be an undifferentiated primitive mantle, consisting of 
about 57% olivine, 17% orthopyroxene, 12% clinopyroxene, and 14% garnet. In this 
model new oceanic lithosphere is formed by the uprise of pyrolite from the low velocity 
layer. It yields a basaltic crust upon partial melting, and leaves a mineralogically 
differentiated and chemically stratified peridotitic lithosphere. The 400- and 650-km 
discontinuity are attributed to phase transformations in the MgO—FeOQ-SiO, system: the 
400-km discontinuity is explained by the divariant phase transformation of olivine to B- 
spinel, and the 670-km discontinuity (originally) by the ’spinel-postspinel’ phase change, 
where the postspinel phase was later experimentally determined as an assemblage of 
perovskite and magnesiowustite. 

The essence of this pyrolite model has not changed, but significantly more detail has 
been suggested concerning the fate of the subducting slab and its role in mantle convection. 
On the basis of the petrology of subducted lithosphere and evidence of phase 
transformations and associated density variations, Ringwood (1982) and Ringwood & 
Irifune (1988) suggest that the harzburgitic and basaltic components of young slabs may 
become trapped on top of the 670-km discontinuity, whereas relatively old slabs penetrate 
into the lower mantle. The 670-km discontinuity is, in this view, a combination of a phase 
transformation with the chemical transition of lenses of buoyantly trapped basaltic crust on 
top of a lower mantle of perovskite and magnesiowiustite. 


PICLOGITE MODEL 


Anderson (1979a, 1984) proposes a chemically stratified upper mantle in which 
eclogite is an important assemblage. He assumes that during the accretional differentiation 
of the Earth heavy eclogite accumulated at the base of the upper mantle, while low-density 
olivine concentrated in the shallow mantle in the form of peridotite. The region of 
accumulated eclogite was originally thought to start at a depth of about 220 km, the 
Lehmann discontinuity (Anderson, 1979a,b), but later (Anderson & Bass, 1984, 1986) it 
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was suggested that this region would coincide with the transition region (between the 400- 
and 670-km discontinuity). Anderson presents several arguments in favour of his model. 
He argues that at an early stage of the history of the Earth melt separation must have 
occured at initially low pressures. From the melt layer garnet crystallized to a 
clinopyroxene-garnet-rich mixture with some olivine, called piclogite (a mixture of eclogite 
and olivine). This assemblage accumulated at the base of the upper mantle, where it is 
gravitationally stable, and thereby created the transition region. The high seismic velocity 
gradient in the transition region is readily explained by the gradual transformation of 
clinopyroxene to garnet-majorite. The 400- and 670-km discontinuity are in this model of 
chemical nature, the 400-km discontinuity being the transition from olivine and 
orthopyroxene to piclogite, and the 670-km discontinuity the transition from piclogite to 
(Mg,Fe)SiO,-perovskite. 


ACCEPTABILITY OF COMPOSITIONAL MANTLE MODELS 


Attempts have been made to ’prove’ the (im)possibility to reconcile these 
compositional models with phenomena of the upper mantle seismic structure. Aspects that 
are of concern in this respect are (1) the general velocity structure in the upper 400 km, in 
the transition zone, and in the lower mantle, (2) the magnitude of the velocity (and density) 
increase at the 400- and 670-km discontinuity, and (3) the width of these discontinuities. 

Although Ringwoods pyrolite model has achieved wide acceptance, there are also 
controversies. Bass & Anderson (1984) and Anderson & Bass (1984, 1986) argued that 
neither the (small) width of the discontinuities, nor the velocity structure within the 
transition zone could be satisfied with the phase transformations in the pyrolitic 
composition. Except for the sharpness of the 670-km discontinuity (originally proposed to 
be the spinel-postspinel phase transformation), these arguments have been attacked by 
Weidner (1985, 1986) and Bina & Wood (1987). In the more recent model (Ringwood & 
Irifune, 1988), the 670-km discontinuity is envisaged as the result of complex interaction of 
subducted lithosphere with the surrounding mantle. The chemical transition of lenses of 
buoyantly trapped basaltic crust on top of the lower mantle could account for the local 
sharpness of the 670-km discontinuity and for the variability of the reflected and converted 
phases as discussed in chapter 5. 

The ’piclogite’ model proposed by Bass & Anderson (1984) and Anderson & Bass 
(1986) has generally been regarded as the counterpart of the essentially uniform pyrolite 
model. The chemical transitions at depths of about 400 and 670 km easily explain the 
sharpness of the discontinuities, although it may be more difficult to explain why converted 
and reflected phases from the 400-km discontinuity are less frequently observed than from 
the 670-km discontinuity (chapter 5). The strong seismic velocity gradient in the transition 
zone is in this model accounted for by the eclogite-majorite transformation. 

Although the laboratory measurements of Knittle et al. (1986) and Knittle & Jeanloz 
(1987) suggest that Mgo 9Fep, ,SiO3-perovskite is more likely as a lower mantle composition 
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(in agreement with the piclogite model) than pyrolite, it is presently not possible to exclude 
either upper mantle model on petrological grounds. A further discussion of the many 
petrological arguments in favour (or against) these compositional models will therefore not 
be given here. 

In conclusion, it should be noted that these upper mantle models should be regarded as 
possible ’simple’ representations of reality: the models are essentially laterally 
homogeneous beneath the low velocity zone (except for regions of subduction or uprising), 
but seismological data, for instance of the European subcontinental mantle, indicate that 
significant small scale lateral heterogeneity is present to a depth of at least 400 km (see 
chapter 4; Spakman et al., 1988). 


6.2 GEODYNAMICS OF THE UPPER MANTLE 
6.2.1 Whole-mantle or layered convection? 


The depth and regime of mantle convection has been a subject of controversy since it 
is accepted that there is mass transport in the mantle. The discussion on this field has 
primarily concentrated on the viability of two extreme models: ’whole-mantle’ or (two-) 
*layered’ convection. In the layered convection model mass transport is confined to two 
separate systems corresponding to the upper and lower mantle, whereas in the whole- 
mantle convection model the system of flow affects the entire mantle. An outline of 
aspects that have contributed to the discussion of whole-mantle or layered convection is 
presented in the following. The observation of a sharp 670-km discontinuity is in this view 
of large relevance, because it indicates that a well-mixed ’whole mantle’ is unlikely. 


INITIAL ARGUMENTS: VISCOSITY OR DENSITY CONTRAST BETWEEN UPPER AND 
LOWER MANTLE? 


Early arguments in favour of ’layered’ (or actually upper mantle) convection were 
based on the assumption of a strong viscosity increase across the transition zone. In initial 
models (e.g. McKenzie, 1967; McConnell, 1968), all mass transport was therefore 
restricted to the upper mantle, while the lower mantle was thought to behave rigidly. 
However, from analysis of postglacial rebound data, it was recognized that there is most 
likely only a modest viscosity increase (at most one order of magnitude) between upper and 
lower mantle (e.g. Cathles, 1975; Peltier, 1976), and this observation has subsequently been 
taken as strong evidence in support of whole-mantle convection (e.g. Davies, 1977; Peltier 
& Jarvis, 1982). 

The notion that the nature of the upper mantle discontinuities may be of relevance to 
possible regimes of mass transport has led several authors to investigate the influence of 
mineralogical phase transformations and compositional density variations on mantle 


The upper mantle 99 


convection. Schubert & Turcotte (1971) and Richter (1973) showed that an exothermic 
phase transition with properties characteristic of the olivine-spinel phase transition 
(associated with the 400-km discontinuity) would not affect the structure of flow, but 
would only increase the flow velocity. An endothermic phase reaction with a Clapeyron 
slope of about -2 MPa/K, as proposed for the spinel-postspinel (670-km) transition (Ito & 
Yamada, 1982), would only have a mildly stabilizing effect (Schubert et al., 1975) and 
would most likely not inhibit whole-mantle convection. This result was later (numerically) 
confirmed by Christensen & Yuen (1985) who estimated that the critical Clapeyron slope to 
force two-layered convection would be in the range of 4 to -8 MPa/K (depending mainly 
on the Rayleigh number of flow). In contrast to the findings for phase transformations, 
Richter & Johnson (1974) showed that if the upper and lower mantle are chemically 
distinct, layered convection would be most probable. 

Another (early) argument for layered convection was therefore based on the inference 
of compositional stratification of the mantle. It was observed that the (seismologically 
determined) density increase across the transition zone was larger than expected from 
phase changes alone (e.g. Anderson, 1968; Press, 1968). However, the empirical 
relationship between density and mean atomic weight on which this argument was based, 
has been found unreliable for high pressure and temperature mineral phases (Ringwood, 
1975). Thus, the question of whether the 670-km discontinuity constitutes a chemical 
transition remained unanswered. 


SLAB PENETRATION INTO THE LOWER MANTLE OR NOT? 


In this respect, evidence from deep focus earthquakes should also be mentioned. 
Studies of seismicity showed that no earthquakes occur below a depth of about 700 km, 
after an increase in seismic energy release between 500 and 700 km (e.g. Isacks & Molnar, 
1971; Vassiliou et al., 1984; Giardini & Woodhouse, 1984). In addition, it was found that 
most deep earthquakes exhibit downdip compressional focal mechanisms (Isacks & 
Molnar, 1971; Giardini & Woodhouse, 1984). These observations have been taken as 
arguments that the lower mantle would act as a barrier to the subducted lithosphere either 
by a strong viscosity increase or a density gradient, which would prevent the slab to 
penetrate into the lower mantle (e.g. Richter, 1979; Vassiliou et al., 1984). The inability of 
slab penetration into the lower mantle would imply two separate mass transport systems in 
upper and lower mantle. 

However, Wortel (1986) showed that the increase in compressional seismic energy 
release between 500 and 700 km and the cessation of seismicity at about 700 km depth can 
also be explained by the thermal assimilation of the subducted plate. He demonstrated that 
no barrier is required to explain the seismicity observations, and this would imply that these 
data may not be invoked as evidence in support of layered convection. 

In fact, there are indications that some subducted slabs do penetrate into the lower 
mantle. Residual sphere analyses (i.e. studies of travel time anomalies that are corrected 
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for all known delays apart from those near the source) of deep earthquakes (Jordan, 1977; 
Creager & Jordan, 1984, 1986; Fischer et al., 1987) and observations of multipathing 
(Silver & Chan, 1986; Beck & Lay, 1986) seem to require aseismic, lower mantle 
extensions of the slabs along the northwestern Pacific of several hundred kilometers. 
Tomographic studies by Grand (1987) and VanderHilst (pers. comm., 1987) show a lower 
mantle S- and P-wave velocity anomaly extending from the eastern Carribean to the 
northern U.S. This anomaly may be associated with past subduction of the Farallon plate 
beneath North America. Although the interpretations of all these studies may be 
questioned, it should be noted that the data are suggestive of mass transport between upper 
and lower mantle in some regions of subduction. 

Thus, if we assume that slab penetration into the lower mantle occurs, does this then 
necessarily imply that the upper and lower mantle are uniform in composition, and is mass 
circulation then ’inevitably’ that of whole-mantle convection? Christensen & Yuen (1984) 
have investigated this aspect by numerical modelling of convection. They showed that, in 
absence of a phase transition, a compositional density contrast of over 5% between upper 
and lower mantle will prevent slabs to descend into the lower mantle. However, if this 
density difference is smaller than 5%, then subducted lithosphere will be able to plunge 
into the lower mantle for several hundred km, and some limited degree of mixing between 
upper and lower mantle will be possible. For an even smaller chemical density contrast, 
less than 2%, the slab will sink to the core-mantle boundary, and the resulting convective 
regime will be essentially that of whole-mantle convection. 


DENSITY CONTRAST BETWEEN UPPER AND LOWER MANTLE 


The density contrast between upper and lower mantle as given by reference Earth 
models is 6% (1066B) to 9% (PREM). As is generally accepted, at least a part of the total 
density increase is due to phase changes. The question is, whether a part of density 
contrast must be contributed to a chemical transition or not. 

Recent experimental work on the thermal expansion coefficient of Mgp Feo ,;SiO,- 
perovskite by Knittle et al. (1986) allowed an estimate of the zero-pressure density of the 
perovskite phase. These authors conclude that the (extrapolated) zero-pressure density of a 
pyrolite at lower mantle conditions is 2% smaller than that estimated for the (seismological) 
lower mantle by adiabatic decompression. Although this does suggest that the upper and 
lower mantle are chemically distinct (assuming a pyrolitic upper mantle), this inference 
may also be questioned, especially concerning to the extrapolation to zero-pressure. In this 
respect, it should be noted that the sharpness of the 670-km discontinuity, implied by the 
observations of P-to-S conversions, is also suggestive of a chemical transition between 
upper and lower mantle. This observations substantiate the (weak) petrological argument 
for a chemical density difference between upper and lower mantle. 

The magnitude of a chemical contribution to the density contrast across the transition 
is not known. If we take, as an example, 2% for the chemical density difference at 670 km, 
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and if we assume that the transition to a predominantly perovskitic phase has a negative 
Clapeyron slope of -2 MPa/K (Ito & Takahashi, 1987), then the modelling results of 
Christensen & Yuen (1984) indicate that slab penetration is possible and that some degree 
of mixing between upper and lower mantle may occur. The style of convection for this 
situation may certainly not be interpreted as whole-mantle convection, but may also not be 
referred to as strict layered convection because of the mass transport between upper and 
lower mantle. Although small variations in the (unknown) compositional density contrast 
or (badly constrained) Clapeyron slope may drastically alter the results, this example 
indicates that the convective regime of the mantle may that of ‘layered convection’ with 
some mass transport between upper and lower mantle. Such a convective regime is 
presently a very acceptable model, because it incorporates phenomena such as a sharp 
670-km discontinuity, a chemical density contrast, and slab penetration. 


GEOCHEMICAL EVIDENCE 


Geochemical contributions to the ’geodynamical discussion’ cannot be excluded in 
this overview, because geochemical data have provided strong arguments in favour of a 
chemically stratified mantle and two-layered convection. 

From analyses of isotope ratios, it was recognized that mid-ocean ridge basalts 
(MORBs) are strongly depleted in incompatible elements, whereas the continental crust 
appeared enriched in these elements. In geochemical models of the mantle this has been 
explained by the ’plate-tectonic process’, which would serve as a mechanism that 
transports incompatible elements from the ridge to the trench, where they are then supplied 
to the continental crust by island-arc volcanism. In this model, the source of MORB and 
the enriched continental crust can be seen as complementary geochemical reservoirs. This 
observation of distinct geochemical reservoirs has led several investigators to estimate the 
size of the depleted reservoir by estimating the mass of the continental crust and assuming 
that the combination of the two must yield an undepleted, primitive mantle composition 
consistent with that of chondrites. Following this approach, Jacobsen & Wasserburg 
(1979), O’Nions et al. (1979), and Turcotte & Kellogg (1986) have reached similar 
conclusions, implying that the depleted reservoir more or less coincides with the upper 
mantle, and that the lower mantle is an undepleted reservoir. If the assumptions are 
correct, this geochemical evidence is a strong argument in support of a chemically stratified 
mantle and two-layered convection. 

There are several uncertainties, however, which may drastically alter the conclusions. 
They concern the way in which early and later differentiation or fractionation has occurred 
in the mantle, the average bulk Earth isotope concentrations, the isotope concentrations and 
volume of the enriched reservoir. Davies (1981) and Allégre et al. (1983) showed that the 
volume of the enriched reservoir may be anywhere between 30 and 90% of the total 
volume of the mantle depending on the assumptions made. Moreover, apart from the 
conventional model’ described above, widely varying convective models have been 
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proposed that may also explain the geochemical observations (e.g. Spohn & Schubert, 
1982; Davies, 1981, 1984; Loper, 1985). 

The observation of a sharp 670-km discontinuity is not consistent with a well-mixed 
uniform mantle. Models with geochemical heterogeneities in the form of distinct ’blobs’ 
moving in a uniform mantle matrix (Davies, 1984) must therefore be excluded. In addition, 
Hoffmann & McKenzie (1985) showed by numerical modelling that only layered 
convection can explain the spatial distribution of geochemical heterogeneities. 

Although it is evident that there are still many controversies regarding the 
interpretation of the geochemical data, Silver et al. (1988) reached a similar conclusion as 
stated in the previous section. They showed that geochemical observations are consistent 
with a model of predominantly layered convection which also includes mass transport 
between upper and lower mantle. 


6.2.2 Geodynamica! models of the upper mantle 


In the previous section, we have concentrated on large scale features of the mantle. 
However, the observations of chapter 4 imply that there is significant small scale lateral 
heterogeneity in the subcontinental mantle beneath Europe. In this section geodynamical 
models of the upper mantle will be discussed that include features of upper mantle 
heterogeneity. 

The key elements of Jordan’s ’tectosphere model’ (Jordan, 1975, 1978) are distinct 
continental roots. These roots consist of mineral assemblages that are less dense than those 
of the suboceanic upper mantle. They serve to stabilize the density variations that arise 
from differences between the oceanic and continental geotherms. These chemical 
variations are thought to be induced by basalt depletion, and this may have caused the so- 
called ’tectosphere’ to extend from less than 100 km (for unstable continents) to 400 km 
depth (beneath stable shields). Seismological data are in support of such a model, as long- 
wavelength regional variations in the seismic velocity structure (extending to a depth of 
300 to 400 km) show a correlation with surface tectonics. The seismic velocity 
distributions of chapter 4 show a high velocity lid that could be interpreted as the 
*tectosphere’ of Jordan’s model. However, our data are also indicative of small scale 
heterogeneity in the upper mantle beneath Europe, and such elements are not incorporated 
in the tectosphere model. 

Viaar (1983) considers the subhorizontal subduction of a young, buoyant oceanic 
lithosphere underneath an old, strong oceanic or continental plate, a mechanism that can 
also be envisaged as the obduction of an old plate over the younger one. In this 
mechanism, ocean ridges are eventually also overridden, a process that presently is 
occuring at the northwestern margin of the U.S., for instance. Vlaar suggests that this 
process of ’lithospheric doubling’ has been important throughout geological history. It 
would induce a strongly layered upper mantle with large (linear) thermal anomalies, and 
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could therefore account for upper mantle thermal anomalies, and for much of the observed 
heterogeneity in magmas. Such a scenario, which is proposed as a mechanism for the 
Alpine orogeny (Vlaar & Cloetingh, 1984), could provide an explanation for the small 
scale lateral heterogeneity inferred from the data presented in chapter 4. However, it is not 
clear to what lateral extent this process may have been important for the upper mantle 
beneath Europe. Although compatible with the heterogeneity predicted with Vlaar’s 
model, the results of chapter 4 lack the horizontal resolution to firmly support or disprove 
the hypothesis of lithospheric doubling. 


6.3 IN CONCLUSION 


Having reviewed the petrological evidence and current theories of mantle dynamics, 
we are able to put the conclusions of chapters 4 and 5 in a somewhat broader perspective. 

Small-scale heterogeneity in the P- and S-wave velocity of the upper mantle beneath 
Europe has been inferred from the analysis of broad-band NARS data. These results 
indicate that the upper mantle is less homogeneous than was generally assumed, and make 
it evident that the classical petrological model for the upper mantle is unsatisfactory. 
Unfortunately, the seismograms from which these results are obtained are too sparse to 
map the three-dimensional P- and S-wave anomalies. A better resolution will be obtained 
from a larger data set that will become available if a much denser network of broad-band 
stations starts operating, as is expected in the next decade. We will then be able to image 
the detailed upper mantle P- and S-wave velocity structure with tomographic techniques 
using waveform information. The research on upper mantle lateral heterogeneity described 
in this thesis can thus be seen as a first step in this direction. Moreover, with the aid of a 
dense network of broad-band stations, we may even be able to map the topography of the 
670-km discontinuity. 

However, further progress in finding more definite answers on questions regarding the 
petrology and geodynamics of the upper mantle should also come from fields outside 
seismology. This is especially pertinent for laboratory investigations of phase 
transformations at conditions such as expected at about 670 km depth. The seismological 
evidence for sharpness of the 670-km discontinuity is now firm, but the absence of hard 
petrological constraints hamper the interpretation. 

The results of this thesis have shown that a detailed analysis of high-quality 
seismological data may contribute new insights into structure of the upper mantle. Small 
scale lateral (seismic velocity) heterogeneity in the European subcontinental upper mantle 
and a (at least locally) sharp 670-km discontinuity are features that must be incorporated in 
petrological or geodynamical representations of the Earth’s upper mantle. 
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Samenvatting 


Seismologisch onderzoek heeft een belangrijke bijdrage geleverd tot onze huidige 
kennis van de structuur van de aarde. De globale verdeling van de seismische snelheden 
was reeds bekend in de 30-er jaren, en op basis hiervan werd een onderverdeling gemaakt 
in verschillende regio’s (korst, bovenmantel, transitie-zone, ondermantel, binnen- en 
buitenkern). De interpretatie van de seismologische gegevens naar de compositionele, 
thermische, en geodynamische opbouw van de aarde is sindsdien een belangrijk onderwerp 
in de geofysica van de vaste aarde. Momenteel staan er nog belangrijke vragen open op dit 
gebied, waarvoor seismologische gegevens met een hoge resolutie van belang kunnen zijn. 
Het onderzoek van dit proefschrift moet in deze context gezien worden, omdat het twee 
aspecten behandelt van de gedetailleerde structuur van de bovenmantel. 

De seismische opbouw van de bovenmantel onder Europa is onderzocht m.b.v. 
ruimtegolven van bevingen uit het oostelijke Middellandse zeegebied die zijn geregistreerd 
in stations van het NARS array in west Europa. Uit analyse van deze gegevens blijkt dat de 
globale seismische structuur van de bovenmantel onder Europa beschreven kan worden 
door een hoge snelheidslaag tot op een diepte van ongeveer 120 km, een lage 
snelheidslaag, en twee discontinuiteiten op dieptes van circa 420 en 650 km. De 
seismische data wijzen echter ook op belangrijke laterale variaties in de snelheidsstructuur 
tot op een diepte van ongeveer 400 km op relatief kleine schaal. Deze resultaten 
impliceren dat de Europese subcontinentale bovenmantel heterogener is dan vaak werd 
verondersteld. 

Een ander onderzoek is verricht naar de gedetailleerde structuur van de 
bovenmanteldiscontinuiteiten. M.b.v. golven die van P naar S zijn geconverteerd in de 
bovenmantel onder de seismische stations, kon worden afgeleid dat de zg. 670-km 
discontinuiteit lokaal bijzonder ’scherp’ moet zijn. Deze observatie kan erop duiden dat de 
boven- en ondermantel chemisch van verschillende samenstelling zijn. 
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